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I.  INTRODUCTION 


High  levels  of  muzzle  blast  overpressure  can  have  adverse  effects  on 
weapon  crew  members,  nearby  structures,  and  Instrumentation,  Control  of  these 
effects  requires  the  ability  to  predict  the  details  of  the  blast  pulse  as  a 
function  of  weapon  design  and  emplacement,  propellant  and  projectile  charac¬ 
teristics,  and  launch  conditions.  Additionally,  techniques  to  reduce  or 
control  the  blast  through  changes  In  these  properties  are  of  practical 
Importance.  Fansler  and  Schmidt^  developed  scaling  relations  that  permit 
estimation  of  peak  Incident  overpressure,  blast  wave  time  of  arrival,  and 
positive  phase  duration.  Comparison  of  these  estimates  with  experiment 
demonstrates  that  the  scaling  approach  provides  a  reasonable  Initial  estimate 
of  muzzle  blast  characteristics.  Since  the  relations  are  sensitive  to  weapon 
launch  conditions  and  propellant  charge  design,  they  may  be  used  to  study  the 
Influence  of  these  properties  on  muzzle  blast. 

The  present  work  has  two  main  objectives.  The  first  Is  to  extend  the 
scaling  approach  to  treat  the  problem  of  pressure  loadings  on  surfaces  adja¬ 
cent  to  the  weapon.  The  second  Is  to  use  the  scaling  relations  In  the  devel¬ 
opment  of  a  computer  code  called  BLAST  that  plots  contour  maps  of  the  muzzle 
blast  quantities.  We  believe  that  the  results  obtained  from  the  sealing 
relations  are  most  easily  Interpreted  when  presented  In  this  form. 

The  remainder  of  this  report  may  be  outlined  as  follows.  First  we  pre¬ 
sent  the  scaling  relations  used  to  calculate  the  Incident  overpressure,  blast 
wave  time  of  arrival,  and  positive  phase  duration  at  points  on  a  surface  that 
can  have  any  desired  orientation  with  respect  to  the  cannon  borellne.  From  the 
calculations  of  Incident  overpressure  and  time  of  arrival,  the  reflected  blast 
overpressure  on  the  surface  of  Interest  Is  determined.  For  all  of  the  blast 
properties,  an  algorithm  Is  established  to  generate  contour  plots.  Represen¬ 
tative  plots  obtained  from  BLAST  are  selected  for  comparison  with  experimental 
results. 


U.  FREE  FIELD  NUZZLE  BLAST  CALCULATIONS 

The  free  field  blast  computations  are  performed  using  scaling  relations 
developed  by  Fansler  and  Schmidt.  Their  scaling  approach  has  been  described 
In  detail;1  hence,  only  the  final  results  of  the  work  are  presented  here. 
The  quantities  of  interest  in  the  Jree  field  blast  prob  at;  ’re  the  peak  inci¬ 
dent  overpressure  in  atmospheres,  P  ,  the  blast  wave  time  cr  rival,  t  ,  and 

a 

the  positive  phase  duration,  t.  The  expressions  derived  for  these  quantities 
are  summarized  by 


V  =  2.4  2 


(la) 


1,  K.  S.  Fansler  and  E.  M.  Schmidt ,  "The  Prediction  of  Gun  Hunsle  Blast 
Properties  Utilising  Sealing  "  U ,  S.  Amy  Ballistic  Reeeai'ch  laboratory t 
Aberdeen  Proving  Ground ,  BEL  Technical  Report  ARBRL-TR-02504 , 

July  1983  (AD  B  0? 58891). 
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where 


t 


a 


(0.94  cose  +  9.24) 


(lb) 


x  -  (i'/aj  [1  +  0.13  (r/i')] 


(lc) 


2  2  1/2 

A'  =  l  [ucose  +  (1  -  y  sin  6)  ] 


(Id) 


Z  -  (r/l') 


-1.1 


,  and 


(le) 


f(Z)  -  1  ♦  1QZ  -  (Z2/1.2)  +  (Z3/2.3)  -  (Z4/3.4)  +  (Z5/4.5)  -  (Z6/5.6).  (If) 


For  subsonic  exit  flow  (  Vp  ^  am  ) 


l 


(8.62  x  10'3)  Pa 
mm 

rv-D 


Cl  ♦ 


1/2 


(2) 


For  supersonic  exit  flow  {  Vp  >  am), 


t  a  (9.28  x  10'2)  D 


Pm  y« 

m  p 


L  (V-lTT-a: 


(1  + 


y(y-1)  V 


2a 


T 


£ 


ro 


(3) 


To  use  these  relations,  It  Is  necessary  to  know  the  weapon 
characteristics  and  propellant  gas  properties  at  shot  ejection.  When  the 
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latter  are  not  known  from  experiment,  they  may  be  determined  using  interior 
ballistics  theory.  The  user  of  BLAST  has  the  option  of  employing  the 
following  simplistic  interior  ballistics  model  to  obtain  the  gas  properties  at 
shot  ejection. 

Lagrangian  interior  ballistics  as  given  by  Corner2  is  utilized.  The 
internal  energy  of  the  propellant  gases  immediately  prior  to  projectile 
ejection  is  given  by 


E  = 


BCRT 


mean 


y=i 


BCRT 

a 

“TT 


-  ltl  (nij  +  C/3)  (1  +  x)  Vp2 


where  B  is  the  fraction  of  propellant  burnt,  C  is  the  propellant  mass,  mj  is 

the  effective  projectile  mass  accounting  for  bore  friction,  and  x  is  the 
ratio  of  heat  losses  to  kinetic  energy.  The  following  semi-empirical 
formulation  for  x  is  used 


(10500)  UD3/2  (Ta-300)  n 

X  3  - • — 9 - TT? - 9 - ffn  (5) 

Ae  Vp^  (m  +  C/3)  [1.7  ♦  671Di/<f  (DVC)*®0] 

where  n  is  the  roughness  factor.  In  this  report  the  roughness  factor  is  taken 
as  1.35,  an  average  of  large  and  small  gun  values.  If  we  assume  an  ideal  gas 
th<*  muzzle  pressure  is  given  by 


P 


m 


(6) 


where  U  is  the  combined  volume  of  the  bore  and  chamber.  The  following 
expression  is  obtained  for  the  exit  sound  speed  of  the  propellant  gases 


affl  -  (yRTm)1/2  -  jy(v-l)  E/[C  ♦  C2/(3ml )]}1/2  *  (7) 


Results  obtained  from  this  simple  theory  were  compared  with  those  of  a  more 
exact  computer  model  3  and  shown  to  produce  similar  results.  Once  the  weapon 


2.  </.  Comer,  Theom  of  the  Interior  Ballierios  of  Gum.  John  Wiley, 
Alev  York,  1950. 

3.  Private  oonnunioation  from  Brian  Bertixind. 
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characteristics  and  launch  conditions  are  determined,  the  scaling  relations 
given  In  Equation  (1)  are  used  to  compute  F,  t,  and  x  as  functions  of  r  and 

a 

e.  The  contour  plotting  algorithm  requires  the  evaluation  of  these  functions 
at  regularly  spaced  points  (gridpoints)  In  the  plotting  domain.  The  values  of 
r  and  9  must  therefore  be  determined  for  each  grldpolnt.  The  geometry  of  the 
problem  Is  shown  In  Figure  1.  The  origin  of  the  polar  coordinate  system  in 
which  contours  are  plotted  Is  the  perpendicular  projection  of  the  muzzle  of 
the  gun  Into  the  plane  of  interest.  The  distance  from  the  origin  to  the 
muzzle  is  h,  and  $  is  the  angle  between  the  plane  and  the  boreline  of  the 
gun.  The  contour  algorithm  requires  a  rectangular  grid  in  the  contour 
plane.  For  this  reason,  the  gridpolnts  are  located  in  a  cartesian  coordinate 
system  originating  at  the  muzzle  of  the  gun.  The  x-y  plane  of  this  coordinate 
system  is  taken  parallel  to  the  contour  plane;  thus,  any  point  In  the  contour 
plane  is  defined  by  (x,y,-h).  The^distance  from  the  muzzle  to  a  point  in  the 
•  plane  is  given  by  the  magnitude  of  r  : 


r 


2  l/Z 
*  IT). 


(8) 


Let  u  be  a  unit  vector  parallel  to  the  boreline.  The  angle  6  between  u  and 
r  is  given  by 


9  =  cos 


-1 


♦ 

r 


+ 

u 


cos 


•1 


ycos$-  h  sln$ 

L(x2  ♦  /  ♦  hZ)I/2 


(9) 


The  values  of  r  and  e  may  be  used  to  obtain  F,  t  ,  and  t  from  the  scaling 

a 

relations  at  each  point  in  the  contour  plane. 


111.  REFLECTED  SHOCK  WAVE  CALCULATIONS 

The  reflection  of  shock  waves  from  surfaces  can  be  quite  complex. 
Edney4  has  identified  a  number  of  possible  flow  structures  depending  upon  the 
strength  of  the  incident  wave  and  the  angle  of  reflection.  In  the  present 
context,  the  process  is  considerably  simplified.  Only  two  types  of  reflection 
are  considered:  regular  and  single  Mach  stem.5 


4.  B.  S.  Rdney ,  Effects  of  Shock  Impifi^emant  on  the  Heat  Tiymofer  Around 
Blunt  Bodies , w A JAA  Jout'mi,  Vol.  16 ,  Bo.  1,  January  1968 1  pp,  16-21, 

5.  B.  P.  Bei'ii'and,  "Measurement  of  Pressure  in  Mach  Reflection  of  Stt'ong 
Shock  Wave3  in  a  Shock  Tube,  "  V.  S.  Amy  Ballistic  Research  Laboratoi'y , 
Aberdeen  Proving  Ground ,  Hai'yland,  BRL-MR-2196t  June  1972  (AD  746613), 
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For  points  very  close  to  the  reflecting  surface  (the  contour  plane),  the 
blast  may  be  considered  planar.  The  shock  wave  angle  of  incidence,  c^,  is 

simply  the  angle  between  the  direction  of  propagation  of  the  blast  and  the 

surface  normal,  n,  to  the  contour  plane  (Figure  2).  One  method  of  obtaining 

the  direction  of  propagation  of  the  blast  wave  is  to  evaluate  the  gradient  of 
ta  since  the  shape  of  the  wave  is  equivalent  to  the  surface  defined  by 

constant  t,.  This  gradient  is  normal  to  the  blast  wave  surface;  therefore,  it 
points  in  the  direction  of  propagation.  The  difficulty  with  this  approach  is 
that  the  computer  time  required  to  evaluate  the  gradient  of  ta  at  each  grid 

point  is  excessive.  In  Appendix  A  we  describe  an  alternate  method  for 
calculating  the  shock  wave  angle  of  incidence.  This  alternate  method  may  be 
shown  to  be  mathematically  equivalent  to  the  gradient  approach,  but  It  Is  more 
efficient  computationally.  In  BLAST  we  use  the  approach  described  in  Appendix 
A. 


Once  the  shock  wave  angle  of  incidence,  e^,  is  obtained,  we  determine 

whether  regular  reflection  or  Mach  reflection  of  the  shock  wave  occurs.  The 
first  step  in  this  procedure  is  to  shift  from  a  fixed  reference  frame  to  one 
moving  along  the  reflecting  surface  at  the  same  velocity  as  the  shock  waves. 
In  this  reference  frame  the  shock  wave  is  stationary  and  has  a  streamline 
flowing  through  it  parallel  to  the  reflecting  surface  (Figure  3).  The 
relations  for  compressible  flow  through  an  oblique  shock  wave6  are  applicable 
to  this  system. 

The  Mach  number  of  the  streamline  in  the  region  in  front  of  the  incident 
shock  is  given  by 


H 


1 


(10) 


where  pj  and  are  the  pressures  behind  and  in  front  of  the  incident  shock, 
respectively.  Across  the  incident  shock,  the  flow  is  deflected  through  an 
angle,  Sj,  given  by 


5 


i 


H  2sin2a.-l 

2(cot  ou  )  -w — - 

Hj  (y*cos  2cij )  ♦  2 


(ID 


S,  H.  V,  Liepeann  and  A .  Roahko,  Elements  of  Gasdunarri cs,  J.  Vi  ley, 
Rev  York,  195?, 
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and  the  Mach  number  behind  the  incident  shock  is  obtained  from 


M 


2 


(Y-l)  (y-1) 

r  00 


1/2 


(12) 


The  question  of  whether  regular  reflection  is  possible  may  now  be 
addressed.  The  maximum  stream  deflection  for  a  specified  upstream  Mach  number 
is  given  approximately  by 


6 

max 


_4 _ 

l/T  (Y+l) 


(13) 


If  2  <$max  »  regular  reflection  is  not  possible  and  Mach  reflection 

occurs.  If  Sj  <  6  ,  regular  reflection  occurs  and  the  analytic  solution  for 

the  reflected  overpressure  may  proceed. 

If  regular  reflection  is  possible,  the  flow  is  deflected  through  the 
reflected  shock.  This  deflection  angle,  5^,  must  be  equal  to  6^  because  the 

boundary  condition  requires  the  flow  to  be  parallel  to  the  reflecting 

surface.  Equation  (11)  applied  to  the  flow  across  the  reflected  shock,  with 

Sj  •«  6?,  yields 


2 

sin 

2  (cot  ou)  — 7 - - -  tan  <5_  *  0  (14) 

c  (y  ^  cos  2^)  c 


where  a9  is  the  wave  angle  of  the  reflected  shock.  Equation  (14)  Is  solved 
for  a->  by  iteration.  Once  is  known,  the  following  equation  is  used  to 
determine  the  pressure  behind  the  reflected  shock,  pq 


pR 


.  2 
sin 


cu-1 )  ♦  1  . 


(15) 


Finally  the  reflected  overpressure  in  atmospheres  is  given  by 


PR  ~  p«  _  £]_ 

P®  P  J  Pco 


(16) 


A  computational  problem  arises  in  the  preceding  development  as 
cij  approaches  zero.  The  expression  for  given  in  Equation  (10)  becomes 

arbitrarily  large  due  to  the  sin  term  in  the  denominator.  To  avoid  this 

difficulty  we  do  not  use  the  relations  presented  above  for  angles  of  incidence 
less  than  one  degree.  For  angles  of  incidence  less  than  this  value,  the 
solution  for  plane  shock  wave  reflection  normal  to  the  surface  is  adequate. 
In  this  special  case  the  ratio  of  reflected  overpressure  to  incident 
overpressure  is  given  by 


pR  (3y-1)  Pj  -  (V"l) 

P7  ~  TV-1)  Pj  +  pit 


If  regular  reflection  does  not  occur  (6.  >  5  )  it  is  impossible  to 

obtain  an  analytic  solution  for  the  reflected  overpressure  using  the  theory  of 
simple  oblique  waves.  In  this  case  the  reflected  overpressure  is  extracted 
from  empirical  results.  A  typical  data  plot  of  reflected  overpressure  versus 
shock  wave  angle  of  incidence  for  various  Incident  overpressures^  is  presented 
in  Figure  4.  Similar  data  from  a  variety  of  sources  were  examined  and  two 
general  features  noted.  First,  beyond  a  certain  angle  of  incidence  whose 
value  depends  on  the  incident  overpressure,  the  reflected  overpressure  curves 
may  be  approximated  by  straight  lines.  To  determine  the  equations  of  these 
lines  a  point  and  slope  are  needed.  At  =  90°,  the  blast  wave  is  at  grazing 

incidence  to  the  surface  and  the  ratio  of  reflected  to  incident  overpressure 
equals  one,  Independent  cf  the  incident  overpressure.  Thus,  the  right  end 
point  of  al*  the  straight  lines  Is  defined.  The  slopes  of  the  lines  are  taken 
to  be  a  function  of  the  incident  overpressure  and  are  obtained  for  11  values 
of  this  parameter,  lagranqiao  interpolation  is  then  used  to  obtain  a  Slope 
for  any  intermediate  value  of  incident  overpressure.  The  Interpolation 
procedure  is  also  used  to  determine  the  angle  of  incidence  beyond  which  the 
linear  approximation  Is  valid. 


7.  S .  Glass  tone,  Editor,  The  Fffectr,  of  Unclear  Weapons.  Dept,  of  the  Army 
Pasrphlet  .Vo.  SO-3,  March  1977,  p.  123. 
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The  second  feature  apparent  in  the  empirical  plots  pertains  to  the  region 
near  the  maxima  of  the  reflected  overpressure  curves.  For  values 

of  cij  between  the  regular  reflection  and  the  linear  interpolation  regimes,  the 

curves  have  a  shape  that  can  be  approximated  by  cubic  polynomials  in  a,.  The 
four  boundary  values  needed  to  evaluate  the  polynomial  are  provided  from  the 
solutions  at  points  A  and  B  of  Figure  5.  The  reflected  overpressure  at  point 

A,  it  (a^),  is  determined  from  the  analytic  solution  for  regular  reflection. 

The  slope  at  A,  ir‘(a^)  ,  can  be  approximated  by  taking  a  finite  difference. 

At  point  B,  Tr(an)  and  tt ‘ (ag)  are  obtained  from  the  Lagrangian  interpolation 
procedure  described  above.  The  four-term  polynomial  at  A  is 


ir^A)  =  %  +  W  W  +  *3°^.  <18) 


The  derivative  of  this  polynomial  is 


»'(aA)  =  *1  +  2lT2aA  +  3lt3aA2  (19) 


with  two  corresponding  equations  at  point  B,  These  four  equations  are  solved 
for  the  four  unknown  coefficients  itg,  tt^ ,  t^,  and  ir^  . 

The  joining  of  the  regular  reflection  solution,  the  cubic  fit,  and  the 
linear  fit  produces  a  continuous  variation  of  reflected  overpressure 
coefficient  versus  angle  of  incidence  for  any  specified  incident 
overpressure.  Results  of  this  procedure  are  plotted  in  Figure  6.  While  this 
method  of  determining  the  reflected  overpressure  is  relatively  simple,  it 
substantially  reproduces  the  empirical  plots  as  is  demonstrated  by  comparison 
with  experimental  data  (Figure  7). 


IV.  THE  CONTOUR  PLOTTING  METHOD 

The  results  :*f  the  analysis  have  been  formulated  as  a  computer  code, 
BLAST,  which  is  programmed  in  FORTRAN  77  and  Is  designed  for  interactive  use 
on  graphics  terminals.  A  listing  of  BLAST  may  be  found  in  Appendix  B. 

The  method  used  to  plot  the  contours  is  relatively  standard  and 
originated  from  an  ALGOL  computer  code.8  A  grid  is  constructed  covering  the 


8.  B.  R.  Heap,  et  at,  "Three  Contouring  Algorithms , "  NPL-81,  National 
Physical  Laboratory Teddington ,  England,  December  1969. 
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region'  of  space  that  Is  of  Interest.  This  grid  is  composed  of  many  small 
rectangular  cells'.  At  each  corner  of  every  cell  (a  gridpoint),  the  value  of 
the  function  to  be  plotted  must  be  known.  Each  contour  value  Is  processed 
individually,  i.e.,  all  contours  at  a  particular  level  are  plotted  before 
moving  to  the  next  level.  The  search  for  contours  at  <;  particular  level 

begins  by  locating  all  the  contours  that  intersect  the  outside  boundary  of  the 
grid  (open  contours).  The  search  along  the  boundary  progresses 

counterclockwise  from  the  bottom  left  corner  of  the  grid  until  two  adjacent 
gridpoints  are'  found  that  bound  the  contour  level.  To  avoid  repeatedly 

finding  and  plotting  the  same  contours,  an  additional  requirement  is 

imposed:  the  larger  functional  values  must  be  to  the  right  of  the 

intersection  point  of  the  crntour  and  the  grid  boundary  for  the  point  to  be 
recognized  as  the  beginning  of  a  contour.  Once  the  beginning  of  a  contour  is 
found,  it  is  followed  cell  by  cell  through  the  grid  by  comparison  of  adjacent 
grid  points.'  The  precise  path  ci  the  contour  across  any  cell  boundary  is 
determined  by  linear  interpolation.  The  coordinates  of  this  interpolated 
point  ahe  passed  to  a  plotting  subroutine  which  uses  a  commercial  graphics 
package,  DISSPLA,  to  draw  each  line  segment  on  the  plotting  device.  Because 
the  grid  squares  are  very  small,  there  is  r:>  perception  of  the  contour  lines 
being  composed  of  straight  segments. 

After  all  the  <■ ar  contours  at  5  particular  level  are  located  and 
plotted,  the  grid  Is  cached  for  closed  contours  at  that  level.  The  closed 
contours,  once  found,  are  followed  through  the  grid  and  plotted  In  the  same 
manner  as  the  open  contours.  When  plotting  of  the  closed  contours  is 
complete,  the  process  is  repeated  for  the  next  contour  level. 


V.  RESULTS 

In  Figures  8  and  9  we  present  contour  maps,  generated  by  BLAST,  of  peak 
incident  overpressure,  peak  reflected  overpressure,  blast  wave  time  of 
arrival,  and  positive  phase  duration.  Maps  of  incident  overpressure,  time  of 
arrival,  and  positive  phase  duration  are  generated  within  thirty  seconds  on  a 
VAX  11/780.  Plots  of  reflected  overpressure  require  additional  computations 
and  are  produced  In  approximately  five  minutes.  The  predictions  shown  are  for 
the  30  mm.  XM230  Chain  Gun.  The  weapon  is  positioned  parallel  to  the  contour 
plane  with  the  muzzle  0.26  meter  distant  from  the  origin  of  the  contour 

grid.  The  strong  directional  dependence  of  the  muzzle  bUst  field  is  readily 
apparent  in  the  plots  of  incident  and  reflected  overpressure.  The  contours  of 
time  of  arrival  represent  the  blast  wave  surface  as  it  propagates  outward  from 

the  muzzle.  This  surface  rapidly  becomes  spherical  with  a  center  displaced  in 

the  forward  direction.  The  positive  phase  duration  map  shows  the  rate  of 

increase  of  the  blast  wavelength  to  be  larger  in  the  forward  direction  than  to 
the  rear  of  the  muzzle. 

To  check  the  validity  of  our  model,  we  compare  BLAST  predictions  with 
experimental  blast  measurements  obtained  in  a  previous  test  program9  of  the  30mm 


9,  h\  M.  Sclvnidt,  "Muzzle  Blast  Pressure  Loadings  on  Aii'craft  Surfaces j  " 
U.  S.  Aimy  Ballistic  Rescai\m  Laboi\iiory ,  Abeixlccn  Pi'oving  Ground , 
Mery  lend t  ARBRL-MR-03336,  Fabi'uamj  1984  (ADA  139132). 
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Chain  Gun.  For  these  tests,  pressure  transducers  were  mounted  In  a  linear 
array  on  an  aluminum  plate.  The  plate  was  then  aligned  at  angles  of  zero  and 
minus  five  degrees  relative  to  the  line  of  fire,  with  separations  of  0.26 
meter  and  0.23  meter  respectively.  For  the  zero  degree  orientation,  data  were 
collected  with  the  transducer  array  aligned  both  parallel  and  normal  to  the 
line  of  fire.  For  the  minus  five  degree  orientation,  only  the  parallel  array 
alignment  was  used.  In  Figures  10,  11,  and  12  we  compare  the  experimental 
measurements  of  reflected  overpressure,  time  of  arrival,  and  positive  phase 
duration  for  each  of  the  three  transducer  orientations.  In  these  plots  the 
abscissa,  S,  denotes  distance  on  the  plate  surface  measured  from  the  point  on 
the  plate  nearest  the  muzzle.  Negative  values  of  S  indicate  positions  behind 
the  muzzle. 

Figure  10  shows  that  the  peak  reflected  overpressure  measurements  agree 
very  favorably  with  BLAST  predictions,  particularly  along  the  line  of  fire. 
The  small  plateau  which  is  apparent  in  each  of  the  predicted  curves  denotes 
the  transition  from  regular  reflection  to  Mach  reflection  of  the  incident 
wave.  For  this  blast  field,  the  effect  appears  too  small  to  be  defined 
experimentally.  From  the  plots  of  time  of  arrival  in  Figure  11,  wo  see  that 
the  location  of  minimum  ta  predicted  by  BLAST  is  shifted  forward  from  the 

location  measured  experimentally;  however,  the  shape  of  the  experimental  ta 

distribution  is  predicted  quite  well.  It  is  obvious  from  Figure  12  that  there 
is  a  definite  deficiency  in  the  positive  phase  duration  scaling  relation.  We 
note  that  the  expression  for  t  was  determined  by  fitting  scaled  experimental 
data  and  that  the  scatter  in  this  data  justified  only  a  simple  linear  curve 
fit.  Clearly,  a  more  complete  analysis  is  desirable. 


VI.  CONCLUSIONS 

The  muzzle  blast  pressure  distribution  on  a  surface  located  In  the 
vicinity  of  a  cannon  has  been  calculated.  The  computational  procedure 
accounts  for  shock  reflection  processes  and  is  based  upon  previously  developed 
scaling  relations  describing  muzzle  blast.  The  approach  Is  formalized  in  a 
computer  program,  BLAST,  that  calculates  and  plots  contour  maps  of  peak 
incident  overpressure,  peak  reflected  overpressure,  blast  wave  time  of 
arrival,  and  positive  phase  duration.  The  contours  are  plotted  In  a  plane 
having  an  arbitrary  orientation  with  respect  to  the  gun  tube.  The 
overpressure  and  time  of  arrival  predictions  compare  favorably  with 
experimental  measurements.  The  positive  phase  duration  is  over-predicted  by 
BLAST  due  to  problems  in  the  scaling  expression  used.  BLAST  Is  completely 
interactive  and  requires  only  minutes  of  computer  (VAX  11/780)  time  to  run. 
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Figure  3.  The  Oblique  Shock  Wave  Solution  for  Reflected  Shock 
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Figure 


INCIDENT  QUERPRESSURE 


*he  Three  Regions  of  a  Typical  Reflected  Overpressure  vs.  Angle  of  Incidence  Curve  (BLAST 
Generated) 


CALCULATIONS  COMPARED  WITH  OTHER  SOURCES 


SHOCK  WAVE  ANGLE  OF  INCIDENCE 

Comparisons  Between  Reflection  Pressure  Results  Obtained  with  BLAST  and  from  Other  Sources 
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Figure  8.  Contour  Maps  of  Peak  Incident  Overpressure  (a)  and  Peak  Reflected 
Overpressure  (b) 
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Figure  10.  Comparison  of  Predicted  Peak  Reflected  Ovorprossure  with 
Experiment 


2 


O  E»G&.  MWL 
PREDICT), 


a.  0°  Plate  Orientation, 
Parallel  to  Line  of 
Fire 


S&CTERS) 


b.  0°  Plate  Orientation, 
Perpendicular  to  Line 
of  Fire 


c,  -5°  Plato  Orientation, 
Parallel  to  Line  of 
Fire 


S$£IER$ 


Figure  11.  Comparison  of  Predicted  Blast  Wave  Time  of  Arrival  with 
Experiment 
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Figure  12.  Coaparison  of  Predicted  Positive  Phase  Duration  with  Exporiaont 
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Analysis  for  Obtaining  the  Shock  Wave  Angle  of  Incidence 

Wt  wish  to  obtain  the  angle  between  the  incident  shock  and  the  contour 

plane.  In  Figure  A-l,  A  is  the  point  of  interest  in  the  contour  plane,  r  is 

the  Rector  directed  from  the  muzzle  to  A,  e  is  the  angle  between  the  boreline 
and  r5  4  is  the  angle  between  the  boreline  and  the  contour  plane,  and  h  is  the 
distance  from  the  muzzle  to  the  contour  plane  along  a  line  normal  to  the 
plane. 

We  define  p  as  the  vector  which  is  directed  from  the  boreline  to  A  and  is 
normal  to  the  shockwave  surface.  We  then  have 


■  t-t 


(A-l) 


where  t  is  a  vector  along  the  boreli  e  as  shown  in  Figure  A-l.  The  angle 
between  r  and  p  is  designated  as  n  .  The  components  of  p  are  found  to  be 


Px=  XA 


(A-2) 


Py  “  yA  -  coos* 


(A-3) 


Pz  =  zA  -  tsin$ 


(A-4) 


The  magnitude  of  is  found  from  the  law  of  sines: 


sin (180°  -  9  -  n)  sin  n 
r  t 


( A-5 ) 


Upon  rearrangement,  Equation  (A-5)  becomes 


r  sin  n 


r  sec  9 


sin(8+n)  tan  0 


tan  n 


-  -  1 


(A-6) 


The  cosine  of  the  angle  of  incidence  Is  given  by 


Vr 


-  .  -  -  +  5  sin  * 

C0S  “1  "  P  *  2, 1/2 


<P«  +  Py  *  P2  > 


(A-7) 


where  n  is  the  unit  normal  to  the  contour  plane.  To  use  Equation  (A-7)  we 
must  first  obtain  £  .  If  n  is  known,  Equation  (A-6)  may  be  employed  to  find 
5  .  The  discussion  below  describes  how  n  may  be  found. 


Referring  again  to  Figure  A-l  we  see  that 


.  1  <dr\ 

tan  "  ■  ‘  r  W, 


(A-8) 


where  (dr/d0L  denotes  the  derivative  of  r  with  respect  to  9  along  a  contour 

of  time  of  arnval,  t,.  To  obtain  (dr/d8)  the  following  identity  is  used: 

a  '•a 


.  at  ata 

-  -H r  >  / 

t  To 

a 


(A-9) 


Recall  that 


£' 


[X  f (z)  +  C1  +  C2  cose] 


(A-10) 


where 


X  n  r/£\  Z  »  X 


-1.1 


and 


£'  =  £*  (0)  . 
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Hence,  the  chain  rule  may  be  applied  to  yield 


Furthermore,  by  differentiation,  we  find 


,3X,  r  dA* 

6  r  (f)2^ 


(«)  .1, 

hr*  V 


6 


ta  dA '  A1  r 
l36  %  “  V  ‘He  "a 

A  00 


) 


sine  , 


and 


■  0 


The  following  result  is  obtained  from  Reference  8: 

dta  t'  (  X1,1  \ 

‘  r.  17770-  j 

Equations  (A-ll)  and  (A-12)  may  now  be  rewritten  as 


(A-ll) 

(A-12) 


(A-13) 

{A— 14) 

(A-15) 


(A— 16 ) 


(A-17) 


and 


3ta 

t,r) 


l  da* 


a<»  d9 


a«>'ta 


4.1 


1  +  X 


—) 

1.1  / 


—  C  sine 
a„  2 


(A- 18) 


,3ta,  1 !  X1,1 

5r  9  °  »4l  +  X1-1 


(A- 19) 
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Substitution  of  Equations  (A-18)  and  (A-19)  into  (A-9)  yields 


Ot  ■  3T  [x  '  ‘a  (1  +  x'ia>:!  +  V(1  +  x"1-1>  s1ne 


where 


at 
”  a 

A' 


Finally  we  obtain 


ta"”“  (A-21) 

In  summary,  to  determine  ,  we  first  compute  n  from  Equation  (A-21). 

Equation  (A-6)  is  then  used  to  obtain  £  ,  which  is  substituted  into  Equations 

(A-2)  through  (A-4)  to  determine  the  components  of  p  .  Finally,  Equation 
(A-7)  Is  used  to  calculate  a^.  This  approach  is  less  time-consuming  than 

direct  computation  of  the  gradient  of  ta  because  we  make  use  of  dta/dX  which 

is  available  to  us  from  previous  analysis. 


Seorsetry  of  the  Shock  Wave  Angle  of  Incidence  Calculation 
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APPENDIX  B  -  LISTING  OF  BLAST 

C  ifiiimiimemHHi 

C  •  BLAST  * 

C 

C - - - - - - - - - - - 

C  BLAST  GENERATES  CONTOUR  MAPS  OF  PEAK  INCIDENT  OVERPRESSURE, 

C  PEAK  REFLECTED  OVERPRESSURE,  BLAST  WAVE  TIKE  OF  ARRIVAL,  AND  POSITIVE 
C  PHASE  DURATION  OF  MUZZLE  BLAST  FROM  CANNON.  CONTOURS  OF  THESE  QUANTI- 
C  TIES  ARE  PLOTTED  IN  A  POLAR  COORDINATE  SYSTEM  WHICH  IS  LOCATED  IN  AN 
C  ARBITRARY  PLANE  WITH  RESPECT  TO  THE  GUN  TUBE.  THE  FREE  FIELD  BLAST 
C  QUANTITIES  ARE  COMPUTED  USING  SCALING  RELATIONS  DEVELOPED  BY  FANSLER 
C  AND  SCHMIDT  ( ABRL-TR-02504) .  THE  REFLECTED  OVERPRESSURE  IS  COMPUTED 
C  USING  THE  ANALYSIS  OF  HEAPS,  FANSLER,  AND  SCHMIDT. 

C  THE  PROGRAM  IS  WRITTEN  IN  FORTRAN  77  WITH  A  FEW  VAX-11  FORTRAN 
C  STATEMENTS.  THE  PLOTTING  IS  DONE  USING  THE  COMMERCIAL  SOFTWARE 
C  PACKAGE,  DISSPLA.  THE  PROGRAM  IS  DESIGNED  TO  BE  RUN  INTERACTIVELY 
C  ON  SEVERAL  TYPES  OF  TERMINALS  WITH  GRAPHICS  CAPABILITY. 

C 

C  THIS  PROGRAM  WAS  DEVELOPED  BY  THE  FLUID  PHYSICS  BRANCH  OF  THE  LAUNCH 
C  AND  FLIGHT  DIVISION  OF  THE  BALLISTIC  RESEARCH  LABORATORY. 

C 

C  CURRENT  VERSION  CREATED  OCTOBER  84 

C - - - 


PROGRAM  BLAST 

PARAMETER(NL=40 , Q=50, P=1 00) 


C - - - 

C  P  AND  Q  ARE  THE  HORIZONTAL  AND  VERTICAL  GRID  DIMENSIONS, 
C  RESPECTIVELY.  NL  IS  THE  MAXIMUM  NUMBER  OF  CONTOURS  THAT 
C  MAY  BE  PLOTTED. 

C - - - 


LOGICAL  INTBALjGRIDUP 

REAL  PSI ANG ( 1 1 ) , ANGLIN (11), PSISLP (11), SLOPES( 1 1 ) 

REAL  TOL.PI.RTOD.DTOR 

REAL  UPBND , LOWBND , YV AR , ATM, GAMMA , AMIN , XYFAC , D , PROPM 
REAL  TRAVEL , BORV , CHAMV , TOTLV , PROJ , PROBB, SPECF , TEMPK, ASUBM 
REAL  PRATIO , AHBSS , VSUBP , GAMMP , MU , ELEANG , PHI , COSPHI , SINPHI 
REAL  RING , USUBY , USYBZ ,X,Y,Z, ZHEIGH ,MAXR 

REAL  R , COSTHE , SINTHE , GM 1 , EXPO , RFACTR , XMASST , DSQ , DONHAF , AREAB 
REAL  DIVFAC , CHIHET , ENEHG , LSCALE , LPRIME , CAPX , CAPZ 
REAL  PB AR ( P , Q) , TSUB A ( P , Q) , TAU ( P , Q) , PREFL ( P , Q) 

REAL  TNUM , TANETA , XCE , PSUBX, PSUBY , PSUBZ 

REAL  ALPHA1 , PRESIN , XI , PPMAX , CONINC , H( 40 ) , PMAX, PMIN , PRMIN 

REAL  PRMAX , TSAMIN , TSAMAX , TAUMIN , T AUMAX , XSTEP , RATIO ( 2 ) 

REAL  TARGX , TARG Y , TARGZ , PBAR 1 ,PBAR2,TSUBA1 ,TSUBA2 
REAL  TAU 1 , TAU2 , TINC , T ( 500 ) , PRESS( 500 ) , TMAX, THIN , TSUBAM 

INTEGER  TRMNUM , BAUD , MAJOPT , NOQRID , MENU , UNIT 
INTEGER  POPT , IOPT , KLOP , KLAP, I , J , ISIGN , ICOUNT , NLEVLS 
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0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 


INTEGER  FINAL, TEST, LIMIT, FOUND ( 40 ) ,ICNT,NPNTS 


CHARACTER  ESC/27/ 

CHARACTER* 1  ANS , ANS 1 , DUMMY , WUNITS 
CHARACTER* 3  ERASE 
CHARACTER* 15  FNAME,WFILE 
CHARACTER»32  GUN 
CHARACTER*38  LABEL1 
CHARACTER*25  LABEL2 
CHARACTER*26  LABEL3 
CHARACTER*79  LABEL6 
CHARACTER*75  LABEL7 
CHARACTER*39  LABEL8 
CHARACTER«49  NAME 
CHARACTER*44  NAME1 
CHARACTER* 45  NAME2 
CHARACTER*44  NAME3 
CHARACTER*42  NAME4 


COMMON/ ONE/ PBAR1 ,TSUBA1 ,TAU1 

C0MM0N/TW0/PBAR2 , TSUBA2 , TAU2 

COMMON/ GRID/ NLEVLS , H , K , OPEN , XM , YM , ISIGN 

COMMON/INC/RINC, MAXR 

COMMON/TRIG/ RTOD , DTOR , PI 

COMMON/ REFL/GAMMA , TOL , UPBND , LOWBND , XI 

COMMON/LAGRAN/PSIANG,PSISLP, ANGLIN, SLOPES 

COMMON/ PRESS/ PBAR , POPT , KLAP 

COMMON/ ATMOS/ ATM 

COMMON  WORK( 12000) ,XPLOT(250) ,YPL0T(250) ,LSAV 


THE  ELEMENTS  OF  ARRAY  ’ANGLIN’  ARE  THE  ANGLES  (RADIANS)  AS  A  FUNC¬ 
TION  OF  INCIDENT  OVERPRESSURE  AT  WHICH  THE  LINEAR  FIT  PART  OF  THE 
REFLECTED  OVERPRESSURE  CURVE  STARTS.  THE  ELEMENTS  OF  ’PSIANG’ 

ARE  THE  INCIDENT  OVERPRESSURES  CORRESPONDING  TO  THESE  ANGLES. 

THE  ELEMENTS  OF  ARRAY  ’SLOPES’  ARE  THE  SLOPES  OF  THE  LINEAR  APPROXI¬ 
MATION.  THE  ELEMENTS  OF  ’PSISLP’  ARE  THE  INCIDENT  OVERPRESSURES 
CORRESPONDING  TO  THESE  SLOPES.  THESE  FOUR  ARRAYS  ARE  USED  IN  THE 
LAGRANGIAN  INTERPOLATION  SUBROUTINE. 


DATA  PSIANO/2. ,5. ,10.,15.,20.,25. ,30., 40. ,50. ,60. ,70./ 

DATA  ANGLIN/1. 29, 1.1 5, 1.06 ,1.01,. 98,. 96,. 93. -91,. 89,. 87,. 86/ 
DATA  PSISLP/2. ,5. ,10. ,15. ,20. ,25. ,30. ,40. ,50. ,60. ,70./ 

DATA  SL0PES/-4. 1 ,-3- ,-2.4,~2. 1 ,-2. ,-2. ,-1 .9,-1 .8,-1 .7,-1 .6, 
-1.5/ 


TOL=. 000001 
P I = B .  1 4 1 5  93 
RTOD=180./PI 
DTOR=PI/ 180. 
UPBND=PI/2 
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0105 

L0WBND=.01 

0106 

YVARsI .013*10. **5 

0107 

ATH= 14. 70 

0108 

GAMMA  =  1.40 

0109 

AMINsO. 

0110 

Cl  =  9.24 

0111 

C2  =  .94 

0112 

GRIDUP  =  .FALSE. 

0113 

0114 

0115 

r 

C  THE  CURRENT  VERSION  HAS  FIVE  TERMINAL  TYPE  DISTINCTIONS. 

0116 

C— — 

0117 

0118 

WHITE{6,») 

0119 

WRITE(6,*)  •  Speoify  the  type  of  terminal  you  are  using.* 

0120 

WRITE(6,«) 

0121 

WRITE (6,*)  *  0  Petrographies . * 

0122 

WRITE(6,*)  *  1  Tektronix  4010  or  4051* 

0123 

WRITE(6,«)  »  2  Tektronix  4014* 

0124 

WRITE(6|*)  '  3  Hewlett  Packard  2623* 

0125 

WRITE(6 .*)  *4  ID  -  100* 

0126 

WRITE(6,») 

0127 

WRITE(6,*)  *  Input  corresponding  integer.' 

0128 

0129 

RE AD*, TRMNUM 

0130 

90 

IF  (TRMNUM  .LT.  0  .OR.  TRMNUM  .GT.  4)  THEN 

0131 

WRITE(6»*)  'Invalid  ohoioe.  Try  again.' 

0132 

READ* .TRMNUM 

0133 

GOTO  90 

0134 

END  IF 

0135 

0136 

IF  (TRMNUM  .EQ.  0  .OR.  TRMNUM  .EQ.  4)  THEN 

0137 

BAUD  a  480 

0138 

ELSE 

0139 

BAUD  =960 

0140 

END  IF 

0141 

0142 

NOGRID =1 

0143 

XYFACal . 

0144 

0145 

105 

CALL  SCLEAR( TRMNUM) 

0146 

WRITE(6 ,*) 

0147 

WRITE(6,«) 

0148 

WRITE(6,*) 

0149 

WRITE(6,«) 

0150 

WRITE(6 ,•) 

0151 

WRITE(6 ,•) 

0152 

WRITE(6 ,*) 

0153 

WRITE(6 ,*)  «  MAIN  MENU  * 

0154 

WRITE(6 , *) 

0155 

WRITE(6.*)  '  1.  Input  or  read  information  for  initial  plots 

0156 

WRITE(6,*)  '  2.  Change  the  input  information.' 
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0157 

WRITE ( 6 , * )  *  3*  Generate  the  plotting  grid  or  read 

the  grid' 

0158 

WRITE(6 ,*)  *  from  a  file.' 

0159 

WRITE(6,*)  '  4.  Plot  contours.' 

0160 

WRITE(6,*)  '  5.  Write  input  information  to  file  for 

future  use 

0161 

WRITE (6,*)  '  6.  Save  the  current  grid  in  a  file  for 

future  use 

0162 

WRITE (6,*)  *  7.  Terminate  the  program.' 

0163 

WRITE(6,*) 

0164 

WRITE(6,*)  'Input  the  corresponding  integer.' 

0165 

READ* ,  MENU 

0166 

0167 

110 

IF  (MENU  .LT.  1  .OR.  MENU  .GT.  7)  THEN 

0168 

WRITE( 6  ,*)  'Invalid  response,  try  again.' 

0169 

READ*,  MENU 

0170 

GOTO  110 

0171 

END  IF 

0172 

0173 

IF  (MENU  .EQ.  1)  GOTO  1000 

0174 

IF  (MENU  .EQ.  2)  GOTO  2000 

0175 

IF  (MENU  .EQ.  3)  GOTO  3000 

0176 

IF  (MENU  .EQ.  4)  GOTO  4000 

0177 

IF  (MENU  .EQ.  5)  GOTO  5000 

0178 

IF  (MENU  .EQ.  6)  GOTO  3500 

0179 

IF  (MENU  .EQ.  7)  GOTO  6000 

0180 

0181 

1000 

CALL  SCLEAR(TRMNUM) 

0182 

0183 

WRITE(6,S) 'You  have  chosen  to  input  information  for  initial  plot 

0184 

WRITE(6 ,*) 

0185 

WRITE (6,*)  'Do  you  wish  to  read  the  information  from' 

0186 

WRITE(6,»)  'a  data  file?  (Y  or  N)» 

0187 

READ(*f*(A1)»)  ANS 

0188 

IF  (ANS  .EQ.  * Y * )  THEN 

0189 

WRITE (6,*)  'Enter  the  name  of  the  data  file  (15  oharaoters)' 

0190 

READ(*f ' (A15) ')  FNAME 

0191 

UNIT  =  1 

0192 

0PEN(UNIT  =  1 ,FILE  =  FNAME , READONLY , STATUS  =  'OLD* 

) 

0193 

ELSE 

0194 

UNIT  =  5 

0195 

END  IF 

0196 

0197 

CALL  SCLEAR(TRMNUM) 

0198 

0199 

WRITE(6,*)  'Enter  gun  nomenclature  (32  oharaoters)' 

0200 

READ (UNIT, '(A32) ')  GUN 

0201 

WRITE(6,«) 

0202 

WRITE(6,*)  'Enter  bore  diameter  of  gun  (meters).' 

0203 

READ(UNIT, *)  D 

0204 

WRITE( 6 , *) 

0205 

WRITE(6,»)  'Do  you  wish  to  use  interior  ballistic' 

0206 

WRITE(6,*)  'formulation?  (Y  or  N)' 

0207 

READ( UNIT , ' ( A1 ) ' )  ANSI 

0208 
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0209 

0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 


IF  (ANSI  .EQ.  »Y')  THEN 
INTBAL  a  .TRUE. 

WRITE(6,*)  'Enter  propellant  mass  (kg).* 

READ (UNIT ,*)  PRCrM 
WRITE(6,«) 

WRITE (6,*)  ' Enter  the  distanoe  the  projectile  * 
WRITE(6,*)  ‘travels  inbore  (meters).* 

REAB(UNIT,»)  TRAVEL 
BORV  =  . 25*PI*D*D*TRAVEL 
WRITE(6,») 

WRITE(6,*)  'Enter  volume  of  chamber  (cubic  meters).' 
READ(UNIT,»)  CHAMV 
TOTLV  =  BORV  +  CHAMV 
WRITE(6,*) 

WRITE(6,*)  ’Enter  projectile  mass  (kg).' 

READ (UNIT,*)  PR0J 
WRITE(6 ,*} 

WRITE(6,*)  'Enter  fraotion  of  propellant  burnt.' 

READ ( UNIT,*}  PR0PB 
WRITE(6,») 

C - - - - - 

C  THE  SPECIFIC  FORCE  IS  EQUAL  TO  THE  GAS  COEFFICIENT  TIMES 
C  THE  ADIABATIC  FLAME  TEMPERATURE  OF  THE  PROPELLANT. 

C - - - 


CALL  SCLEAR(TRMNUM) 

WRITE ( 6, •)  'Enter  specific  force.  Use  MKS  units  (m*m/s*s)  (A),' 
WRITE(6,*} 'or  English  units  (ft*lb/lbm)  (B).  Enter  A  or  B.' 
WRITE(6 ,*) 

READ (UNIT | ' ( A1 ) » )  WUNITS 
WRITE(6 ,*) 

WRITE(6,*)  'Now  enter  specific  foroe.’ 

WRITE (6,*) 

READ(UNIT, *)  SPECF 

IF  (WUNITS  .EQ.  »B»)  SPECF  =  SPECF*2.992 
WRITE (6,*) 

WRITE (6,*)  'Enter  flame  temperature  (deg.  Kelvin).' 

’  WRITE(6,*)  'Default  =  3000  K.» 

READ(UNIT,'(F12.4)»)  TEMPK 
IF  (TEMPK  .EQ.  0.0)  TEMPK  a  3000.0 
CALL  SCLEAR(TRMNUM) 

ELSE'  " 

INTBAL  =  .FALSE. 

WRITE(6,*) 

WRITE (6,*)  'Enter  pre-unoorked  sound  speed  at  muzzle  (m/s).' 

READ  (UNIT,  *>)  ASUBM 

WRITE(6,*) 

WRITE(6,*)  'Enter  ratio  of  pre-unoorked  muzzle  pressure' 
WRITE(6,*)  'to  ambient  pressure.' 

READ (UNIT,*)  PRAYIO 
END  IF 
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0261 

0262 

0263 

0264 

0265 

"0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 

0307 

0308 

0309 

0310 

0311 

0312 


WRITE(6,») 

WRITE(6,*)  'Enter  ambient  sound  speed  (m/s).  Default  a  340.0  m/a' 
READ (UNIT, '(F12.4)')  AMBSS 
IF  (AMBSS  .EQ.  0.0)  AMBSS  a  340.0 

WRITE(6,») 

WRITE(6 ,*)  'Enter  exit  velooity  of  projeotile.' 

READ (UNIT,*)  VSUBP 

WRITE( 6 , * ) 

WRITE(6,*)  'Enter  propellent  speoifio  heat  ratio,  gamma.' 

WRITE(6 ,*)  'Default  =1.24.' 

READ ( UNIT , ' ( F7 • 4 ) ' )  GAMMP 
IF  (GAMMP  .EQ.  0.0)  GAMMP  a  1.24 

WRITE(6 ,*) 

WRITE(6,*)  'Enter  momentum  index,  mu.  Default  a  0.78' 

READ ( UNIT , » ( F7 . 4 ) * )  MU 
IF  (MU  .EQ.  0.0)  MU  =  0.78 


CALL  SCLEAR(TRMNUM) 

WRITE(6 ,*)  'Enter  elevation  angle  of  gun  (degrees).' 
READ (UNIT,*)  ELEANQ 
PHI  =  DT0R*ELEANG 
COSPHI  =  COS (PHI) 

SINPHI  a  SIN (PHI) 

WRITE(6 ,*)  'Enter  the  height  of  the  muzzle  (meters).' 
READ (UNIT,*)  ZHEIGH 

CLOSE  (UNITal) 


GOTO  105 

2000  CALL  SCLEAR(TRMNUM) 


'You  have  ohosen  to  change  information.' 
'  Change  of  information  menu' 


2010 


WRITE(6 , *) 

WRITE(6 ,*) 

WRITE(6 ,*) 

WRITE(6 ,*) 

WRITE(6,*) 

WRITE(6 ,*) 

WRITE(6 ,*) 

WRITE(6 ,*) 

WRITE(6,«) 

READ*, MENU 

IF  (  MENU  .LT.  1  .OR.  MENU  .G7.  3  )  THEN 
WRITE(6,*)  'Invalid  response.  Try  again.' 
READ*, MENU 
GOTO  2010 


1 . 
2. 
3. 


Change  weapon  and  flow  parameters.' 
Change  the  orientation  of  the  gun.' 
Return  to  the  main  menu.' 


'Enter  the  corresponding  integer.' 
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0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 

0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 

0343 

0344 

0345 

0346 

0347 

0348 

0349 

0350 

0351 

0352 

0353 

0354 

0355 

0356 

0357 

0358 

0359 

0360 

0361 

0362 

0363 

0364 


IF  (MENU  .EQ.  1)  GOTO  2100 
IF  (MENU  .EQ.  2)  GOTO  2200 
IF  (MENU  .EQ.  3)'  GOTO  105 

2100  CALL  SCLEAR(TRMNUM) 

WRITE ( 6, a)  'You  may  change  any  of  the  following:' 

WRITE(6,») 

WRITE(6,a)  '  1.  Gun  nomenclature  10.  Pre-unoorked  sound  speed 
WRITE (6,*)  *  2.  Bore  diameter  at  the  muzzle  ' 

WRITE(6,*)  '  3.  Propellant  mass  *  11.  Ratio  of  pre-unoorked  muz 

WRITE (6,*)  *  4.  Inbore  travel  zle  pressure  to  ambient 

WRITE(6,*)  •  distance  *  pressure  ' 

WRITE (6,*)  '  5.  Chamber  volume  •  12.  Ambient  sound  speed' 

WRITE(6,*)  '  6.  Projectile  mass  *  13.  Projectile  exit  velocity' 

WRITE (6,*)  •  7.  Fraction  of  pro-  14.  Speoifio  heat  ratio' 

WRITE(6,*)  '  pellant  burnt  •  15.  Momentum  index* 

WRITE (6,*)  '  8.  Speoifio  foroe  •' 

WRITE (6,*)  '  9.  Flame  temperature  •* 

WRITE(6,*) 

WRITE(6,*)  '  or  you  may  16.  Return  to  proceeding  menu 

WRITE(6 ,*) 

WRITE(6,*)  '•  indicates  an  interior  ballistics  parameter.* 
WRITER,*} 

WRITE(6,a)  'Enter  corresponding  integer.* 

READ*, MENU 

2110  IF  (MENU  .LT.  1  .OR.  MENU  .GT.  16)  THEN 

WRITE£6,*)  'Invalid  response.  Try  again.* 

READ* , MENU 
GOTO  2110 
END  IF 

CALL  SCLBAR(TRHNUN) 

IF  (MENU  ,EQ.  1)  THEN 

WRITE(6 ,*)  'Current  value  is  ',GUN 
WRITE (6,*)  'Enter  new  value* 

READ(*,'(A32)')  GUN 
END  IF 

I**  (MENU  .EQ.  2)  THEN 

WRITE(6,*)  'Current  value  i3  *,D 
WRITE(6,*)  'Ente**  new  value* 

READ*,D 
END  IF 

IF  (MENU  .EQ.  3)  THEN 

WRITE(6 ,*)  'Current  value  is  '.PR0PM 
WRITE (6,*)  'Enter  new  value* 


0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

037*1 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

038*1 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

039*1 

0395 

0396 

0397 

0398 

0399 

0*100 

0*101 

0*102 

0*103 

0*10*1 

0*105 

0*106 

0407 

0408 

0409 

0410 

0411 

0412 

0413 

0414 

0415 

0416 


READ* , PROPM 
ANSI  =  'Y' 

END  IF 

IF  (MENU  .EQ.  4)  THEN 

WRITE (6,*)  'Current  value  is  \ TRAVEL 
WRITE (6,*)  'Enter  new  value' 

READ*, TRAVEL 

BORV  =  0.25*PI*D«D*TRAVEL 
ANSI  =  'Y ' 

END  IF 

IF  (MENU  .EQ.  5)  THEN 

WRITE (6,*)  'Current  value  is  ' ,CHAMV 
WRITE (6,*)  'Enter  new  value* 

READ*, CHAMV 
ANSI  =  *Y' 

END  IF 

IF  (MENU  .EQ.  6)  THEN 

WRITE (6,*)  'Current  value  is  *,PROJ 
WRITE(6,*)  'Enter  new  value* 

READ* , PROJ 
ANSI  =  'Y ' 

end  if 

IF  (MENU  .EQ.  7)  THEN 

WRITE(6 ,*)  'Current  value  is  \PR0PB 
W8ITE(6,*)  ’Enter  new  value' 

READ*,PROPB 
ANSI  «  'Y' 

END  IF 

IF  (MENU  .EQ.  8)  THEN 

WRITE(6 ,e )  'Current  value  is  ’ ,SPECF 

WRITE(6 ,*)  'Use  KKS  units  (m*»/s*s)  (A),  or  English' 

WRITE(6,*) 'units  (ft»lb/lbis)  (B).  Eater  A  or  B.' 

WRITE(6,*) 

READ(*, ' ( A1 ) ' )  WUNITS 
WRITE(6,*) 

WRITE(6,*)  'Now  enter  new  value.' 

READ*,  SPECF 

IF  (WUNITS  . BQ.  1 B' >  SPECF  a  SPECP*2.992 
ANSI  a  'Y' 

END  IF 

IF  (MENU  .EQ.  9)  THEN 

WRITE(6,*)  'Current  value  is  *,TEHPK 
WRITE (6,*)  'Enter  new  value' 

READ*,TEMPK 
ANSI  =  'Y' 

END  IF 
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[IE 


0418 

0419 

0420 

0421 

0422 

0423 

0424 

0425 

0426 

0427 

0428 

0429 

0430 

0431 

0432 

0433 

0434 

0435 

0436 

0437 

0438 

0439 

0440 

0441 

044? 

0443 

0444 

0445 

0446 

0447 

0448 

0449 

0450 

0451 

0452 

0453 

0454 

0455 

0456 

0457 

0458 

0459 

046C 

0461 

0462 

0463 

0464 

0465 

0466 

0467 

0468 


IF  (MENU  .EQ.  10)  THEN 

WRITE (6,*)  'Current  value  is  ',ASUBM 
WRITE (6,*)  ’Enter  new  value’ 

READ® , ASUBM 
END  IF 


IF  (MENU  «EQ.  11)  THEN 

WRITER*)  ’Currant  value  is  *,PRATIO 
WRlTE(6,*)  'Enter  new  value’ 
READ*,PRATI0 
END  IF 


IF  (MENU  -EQ.  12)  THEN 

WRITE(6 ,*)  'Current  value  is  ’,AMESS 
•WRITE(6j*)  ’Enter  new  value' 

1  ..  READ*,AMBSS 
'  -  END  IF 


IF  (MENU  '.EQ.  13)  THEN 
'  .WRITE ( 6 ■, * )  ’Current  value  is  *  ,VSUBP 
WRITE(6 ,*)  'Enter  new  value’ 
READ»,VSUBP 
END  IF 


IF  (MENU  .EQ.  14)  THEN 

WI;ITE(6,*)  'Current  value  is  ',GAMMF 
WRITE (6,*)  'Enter  new  value* 

HEAD* , GAMMP 
END  IF 


IF  (MENU  .EQ.  15)  THEN 

WRITE(6 ,*)  'Current  value  is  ',MU 
WRITE (6,*)  'Enter  now  value' 
READ*, MU 
END  IF 


IF  (MENU  .EQ.  16)  GOTO  2000 
GOTO  2100 


2200  CALL  SCLEAR(TRMNUM) 

WRITE(6,*)  'Current  elevation  angle  of  gun  ia  ’,ELEANG, 
4  '  degrees.' 

WRITE(6 ,*) 

WRITE(6 ,*)  'Enter  the  new  elevation  angle.' 

READ*,  ELEANG 
PHI  a  DTOR*KLEANG 
COSPHI  a  COS (PHI) 

SINPHI  =  SIN (PHI) 


WRITE(6,*) 

WRITE(6,*)  'The  current  height  of  the  muzzle  ic  ',ZHEIGH 


WRITE(6,®) 

WRITE(6,®)  'Enter  the  new  height  of  the  muzzle, ' 

READ®,  ZHEIQH 

GOTO  2000 

C - - - - - - - - - 

C  IN  SECTION  3000  THE  USER  CHOOSES  WHICH  QUANTITY  HE  WANTS  PLOTTED. 

C  THE  VALUES  OF  THIS  QUANTITY  ARE  THEN  DETERMINED  AT  EACH  GRID  POINT 
C  IN  THE  CONTOUR  PLANE. 

3000  CALL  SCLEAR(TRMNUM) 


C  THE  USER  MAY  GENERATE  A  GRID  TO  PLOT  ALL  QUANTITIES,  OR  HE  MAY 
C  CHOOSE  TO  EXCLUDE  THE  REFLECTED  OVERPRESSURE,  THIS  CHOICE  IS 
C  GIVEN  BECAUSE  THE  REFLECTED  OVERPRESSURE  GRID  IS  MUCH  MORE  TIME- 
C  CONSUMING  TO  GENERATE  THAN  THE  OTHER  QUANTITIES. 

C  THE  USER  MAY  ALSO  READ  A  GRID  THAT  WAS  PREVIOUSLY  GENERATED. 


WRITE (6 
WRITE (6 
WRITE (6 
WRITE (6 
WRITE (6 
WRITE (6 
WRITE(6 
WRITE(6 
WRITE (6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE(6 
WRITE  (6 
WRITEC6 
WRITE(6 
WRITEC6 
READ  ®, 


,*) 

,*)  'You  may' 
,•) 


1 .  Generate  a  grid  of  values  for  ' 

psak  inoident  overpressure,' 
peak  refleoted  overpressure, ' 
blast  wave  time  of  arrival,' 
positive  phase  duration,' 

2.  Generate  a  grid  of  values  for  ell  of  the' 
quantities  above  exoluding  peak  refleoted' 
overpressure. ' 

3.  Read  a  previously  generated  grid  fro#  a' 
data  file. ' 

Enter  the  corresponding  integer.  ' 


IF  (IOPT  .HE.  3)  THEN 


C  NEW  PAGE. 
C - 


CALL  SCLEAR(TRHNUH) 

KLOPaO 

KLAPsO 


0521  C  THE  OSES  DEFINES  THE  MAXIMUM  RADIUS  TO  BE  PLOTTED. 

0522  C - - - * - — - - - — 

0523 

0524  .  WRITE(6,*) 'Input  maximum  distanoe  from  origin  to  be  plotted 

0525  A  (meters).’ 

0526  ..WRITER,?) ’The  origin  is  the  intersection  point  of  the  contour 

052?  &  plane  and' 

0528  WRITE(6,*) ’the  line  that  passes  through  the  muzzle  of  the  gun 

0529  4  and’ 

0530  WRITER,*)  ’is  perpendicular  to  the  contour  plane.’ 

0531  WRITE(6,») 

0532  READ  *,MAXR 

0533 

0534  C - — — - - - - - - - - - - - - - 

0535  C  PHI  IS  THE  ANGLE  BETWEEN  THE  BORELINE  AND  THE  PLANE  IN  WHICH 
0536  C  THE  CONTOURS  ARE  TO  BE  PLOTTED.  USUBY  AND  USUBZ  ARE  THE  Y  AND  Z 

0537  C  COMPONENTS  OF  THE  UNIT  VECTOR,  U,  PARALLEL  TO  THE  BORELINE  OF  THE 

0538  C  GUN.  THE  FUNCTION  TO  BE  PLOTTED  IS  EVALUATED  AT  EACH  GRIDPOINT. 

0539  C  THIS  PROCESS  PROCEEDS  ROW  BY  ROW  STARTING  IN  THE  LOWER  LEFT  CORNER 

0540  C  OF  THE  GRID.  THE  GRID  IS  A  RECTANGLE  WHICH  EXTENDS  FROM  0  TO  TWO 
0541  C  TIMES  MAXR  HORIZONTALLY  AND  FROM  0  TO  MAXR  VERTICALLY . THE  MAXIMUM 

0542  C  RADIUS  .IS -ROUNDED  TO  HUNDREDS  IN  ORDER  TO  FIT  ON  THE  PLOT  PROPERLY. 

0543  C - - - - - 

0544 

0545  MAXR=INT(MAXR*P)/P 

0546 

0547  WRITE(6 , *) 

0548  WRITE (6,*) 

0549  WRITE(6 , •) 

0550  WRITER,*) 

0551  WRITE(6 ,*) 

0552  WRITE (6 ,*) 

0553  WRITE(6, *) 

0554  IP  (IOPT.EQ. 1 )  THEN 

0555  WRITE(6,*) ’Working.  The  calculation  of  functional  values  will 

0556  4  take  several  minutes. ' 

0557  WRITE(6 ,*) 'Wait  for  bell.' 

0558  END  IF 

0559  IF  (IOPT  .EQ.  2)  THEN 

0560  WRITE(6,*) ’Working.  Wait  for  bell.' 

056 1  END  IF 

0562 

0563  RINC=2.*HAXR/P 

0564 

0565  USUBY=COSPHI 

0566  USU3Z=SINPHI 

0567 

0568  DO  3160  J=1 ,Q 

0569  DO  3165  1=1, P 

0570 

057 1  C - - - 

0572  C  X,  Y,  AND  7.  ARE  THE  X,  Y  AND  Z  COMPONENTS  OF  THE  VECTOR,  R, 
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i 


*{ 

y 


,1 


0573 

0574 

0575 

0576 

0577 

0578 

0579 

0580 

0581 

0582 

0583 

0584 

0585 

0586 

0587 

0588 

0589 

0590 

0591 

0592 

0593 

0594 

0595 

0596 

0597 

0598 

0599 

0600 

0601 

0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 

0611 

0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 

0621 

0622 

0623 

0624 


C  WHICH  POINTS  FROM  THE  ORIGIN  OF  THE  CONTOUR  PLANE  TO  THE  POINT  IN 
C  THE  CONTOUR  GRID  UNDER  CONSIDERATION.  X,Y,  AND  Z  ARE  THE  COORDI- 
C  NATES  OF  THIS  GRIDPOINT  REFERRED  TO  AN  ORIGIN  AT  THE  MUZZLE.  THETA 
C  IS  THE  POLAR  ANGLE  OF  THE  VECTOR  R.  THETA  IS  DETERMINED  BY 
C  TAKING  THE  DOT  PRODUCT  OF  R  AND  U. 

C - - - - - - - - - - - 


X=-J*RINC 

IF  (I»RINC.LT.MAXR)  Y=-(MAXR-I"RINC) 
IF  (I»RINC.EQ.MAXR)  YcO. 

IF  ( I*RINC . GT . MAXR )  Y=I»RINC-MAXR 
Zs-ZHEIGH 

R=SQRT(X**2  +  Y**2  +  Z»«2) 

COSTHE= ( ( Y*USUBY ) + ( Z*USUBZ ) ) /R 
SINTHE=SQRT( 1 .  -  C0STHE**2) 
GM1=GAMMP-1 . 

EXPO= ( 3 • »GAMMP- 1 . ) /GM1 


C  THE  SCALING  LENGTHS,  LSCALE  AND  LPRIME,  ARE  DETERMINED, 
C  FIRST,  INTERIOR  BALLISTIC  PARAMETERS  ARE  CONSIDERED, 

C- - - - - - - , - , - - 


4 

& 

& 


IF  (INTBAL)  THEN 
PROJM  =  1.05»PR0J 
RFACTR  =1.35 
IF(D.LT,0.02)  RFACTR  =  1.5 
IF(D.OT. 0.045)  RFACTR  =  1.3 
XMASST  =  PROJM  +  PR0PM/3. 

DSQ  =  D»"2 
DONHAF  =  SQRT(D) 

AREAB  =  PI*(D""2)/4. 

DIVFAC  a  1.7  +  671 J#nONHAF*((DSQ/PROPM)#i.86) 

CHIHET=1 0500 . *TOTLV*D»DONHAF»(TEMPK-30O . ) "RFACTR/ AREAB 
/(VSUBP**2)/XMASST/DIVFAC 
ENERQ=(  ( PR0PB"PR0PM»SPECF)/QH1 )  -  .5*(  PROJM*. 

PR0PM/3 . )*( 1 .0+CHIHET)"VSUBP«"2 

ASUBM=  ( { GAMMP*GM  1  "ENERG )  /  ( PROPMt-  ( PR0PM*"2/  ( 3 .  "PROJM 

))))**. 5 

PRATIOs ( GM 1 "ENERG ) / ( TOTLV" ( 1+PROPM/ ( 3 . "PROJM) ) )/YVAR 
END  IF 


C  NEXT  MUZZLE  EXIT  CONDITIONS  FOR  SUBSONIC  EXIT  FLOW. 
C - — - - 


& 

A 


IF  ( VSUBP . LE . ASUBM )  THEN 

LSCALE=D*SQRT( ( .00862"PRATIO"ASUBH/(GM1"AHBSS) ) 
•(  1 . +0AMMP"QM1/  2 . )  •(  ( 2 .  ^*0M1  "VSUBP/ ASUBH)  / 

( QAHMP+ 1 . ) ) ""EXPO ) 

ELSE 
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0625 

0626 

0627 

0628 

0629 

0630 

0631 

0632 

0633 

0631* 

0635 

0636 

0637 

0638 

0639 

0640 

0641 

0642 

0643 

0644 

0645 

0646 

0647 

0648 

0649 

0650 

0651 

0652 

0653 

0654 

0655 

0656 

0657 

0658 

0659 

0660 

0661 

0662 

0663 

0664 

0665 

0666 

0667 

0668 

0669 

0670 

0671 

0672 

0673 

0674 

0675 

0676 


C - — - 

C  FINALLY,,  MpZ'ZLE  EXIT  CONDITIONS  FOR  SUPERSONIC  EXIT  FLOW. 
C - - - - - - - .... _ 


.  LSCALEsD«.0928«SQRT( (PRATI0»VSUBP/(GM1«AMBSS) )* 
Cl • +GAMMP*GM1 * . 5» ( VSUBP/ ASUBH) *»2 ) ) 

END  IF 


C- _ _ _ 

C  LPRIME  IS  CALCULATED. 

C- — - - - - 

LPRIHE=LSCALE* ( (MU«C0STHE)+SQRT( 1 .-(MU**2)«(SINTHE»«2) ) ) 

CAPX=R/LPRIME 

CAPZ=CAPX**-1 . 1 

IF  (CAPX.LT. 1 .)  KL0P=1 


C - - - - - - 

C  PBAR  IS  THE  PEAK  INCIDENT  OVERPRESSURE,  TSUBA  IS  THE  BLAST  WAVE 
C  TIME  OF  ARRIVAL,  AND  TAU  IS  THE  POSITIVE  PHASE  DURATION.  THESE 
C  QUANTITIES  ARE  CALCULATED  USING  THE  SCALING  RELATIONS  DERIVED  BY 
C  FANSLER  AND  SCHMIDT. NEW  PBAR  AND  TSUBA  EQUATIONS  ADDED  MARCH, 83. 
C - 


■  PBAR(I, J)=2.4*CAPZ 

TSUBA(I,J)=R/AMBSS»(1.+10.»CAPZ~.8333#CAPZ«2+0.4348» 
&  CAPZ»»3- .2941 •CAPZ#i4+ . 2273*CAPZ##5- . 1 667*CAPZ»*6 ) - 

4  ( LPRIME/ AMBSS ) * ( Cl +C2*C0STHE ) 

TAU ( I, J) =( LPRIME/ AMBSS)*( 1 .0  +  .13«CAPX) 


C - - - - - - - - - 

C  IF  IOPTsI  WE  GO  THROUGH  THE  REFLECTED  OVERPRESSURE  CALCULATION. 

C  THE  FIRST  STEP  IS  TO  DETERMINE  THE  SHOCK  WAVE  ANGLE  OF  INCIDENCE 
C  AT  EACH  GRID  POINT. 


IF  (IOPT  .EQ.  1)  THEN 


C - - - - - 

C  ALPHA 1  IS  THE  SHOCK  WAVE  ANGLE  OF  INCIDENCE  AT  THE  GRIDPOINT. 

C  A  DETAILED  EXPLANATION  OF  THE  ALPHA 1  CALCULATION  IS  IN  APPENDIX 
C  A  OF  THE  REPORT. 


TSABAR  =  TSUBA( I , J ) •AMBSS/LPRIHE 

C - - - ' - ... 

C  DLPDTH  a  DERIVATIVE  OF  LPRIME  WITH  RESPECT  TO  THETA 
C  DCXDTB  =  DERIVATIVE  OF  CAPX  WITH  RESPECT  TO  TSUBA  BAR 
C - 
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uo77 

0678 

3679 

0680 

0681 

0682 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 
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0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 

0711 

0712 

0713 

0714 

0715 

0716 

0717 

0718 

0719 

0720 

0721 

0722 

0723 

0724 

0725 

0726 

0727 

0728 


DLPDTH  =  -LPRIME*MU*SINTHE/SQRT ( 1 .  -  MU»MU»SINTHB*«2) 
DCXDLP  =  -R/LPRIME»»2 


C - - - — — - - - 

C  IF  THE  TIME  OF  ARRIVAL  EXPRESSION  IS  CHANGED,  DTBDCX  MUST  ALSO  BE 
C  CHANGED  ACCORDINGLY. 

DCXDTB  s  1.  +  CAPZ 

TANETA  =  (DLPDTH/LPRIME) *(  DCXDTB»TSABAR/CAPX  -  1.0  ) 

4  +  C2#LPRIME*DCXDTB*SINTHE/CAPX 

IF  (COSTHE  .EQ.  0.)  THEN 
XCE  =  R*TANETA 
ELSE 

TANTHE  =  SINTHE/COSTHE 

XCE  s  { R/ COSTHE ) / ( TANTHE/TANETA  ♦  1.0) 

END  IF 


PSUBX  =  X 

PSUBY  s  Y  -  XCE*COSPHiI 
PSUBZ  =  Z  -  XCE*SINPHI 

ALPHA  1  =  ACOS  ( -PSUBZ/ SQRT  ( PSUBXM2+PSUBYa*2+PSUB2**2)  ) 
IF  (ALPHA1  .GT.  PI/2,)  ALPHA 1  a  PI  -  ALPHA) 
PRESIN=PBAR(I,J) 

XIaPRESIN+1 . 


C- 

C 

C 

C 

C 

C 

C- 


SUBROUTINE  REFCAL  IS  CALLED  TO  DETERMINE  THE  REFLECTED 
OVERPRESSURE,  PREFL ,  AT  THE  GRIDPOINT. 

IF  THE  ANGLE  OF  INCIDENCE  IS  GREATER  THAN  1 .  DEGREE  HE  USE 
OBLIQUE  SHOCK  THEORY.  IF  ANGLE  OF  INCIDENCE  IS  LESS  THAN 
1.  DEGREE  HE  USE  NORMAL  SHOCK  RELATIONS. 


IF  ( ALPHA 1  .GT.  0.0174)  THEN 

CALL  REFCAL ( PRESIN ,PRBFL( I, J),ALPHA1) 

ELSE 

P3  =  (  XI»«2»(3*GAHHA-1.)  -  {GAKMA-1 .)»XI  )/ 
4  (  XI#(QAHHA-1 . )  ♦  (GAHMA+1.)  ) 

PREFLU, J)  n  P3  -  1. 

END  IF 


END  IF 

3165  CONTINUE 
3160  CONTINUE 
CALL  BELL 
CALL  BELL 
GRIDUP  =  .TRUE. 
GOTO  105 
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0732 

0733 

0734 

0735 

0736 

0737 

C  IF  THE  USES  HAS  CHOSEN  TO  READ  A  GRID  THAT  HAS  BEEN  SAVED 

C  THE  QHID  IS  NOW  READ  FROM  THE  DATA  FILE. 
r. _  _ 

3200 

CALL  SCLEAR(TRMNUM) 

0738 

WRITE (6,*)  'What  is  the  name  of  the  data  file?' 

0739 

READ ( • , '  ( A20 ) » )  FNAME 

0740 

OPEN (UNlTs 1 , FILEsFNAME , STATUS* ' OLD • ) 

0741 

WRITE(6f *)  'Reading  grid  from  file  », FNAME 

0742 

READ( 1 ,*)  IOPT 

0743 

READ( 1  ,•)  MAXR 

0744 

IF  (IOPT  .EQ.  1)  THEN 

0745 

DO  3210  1=1, P 

0746 

DO  3215  J=1,Q 

0747 

READ( 1 ,»)  PBAR(I, J) ,PREFL(I, J) ,TSU3A(I, J) ,TAU(I, J) 

0748 

3215 

CONTINUE 

0749 

3210 

CONTINUE 

0750 

END  IF 

0751 

.  , 

IF  (IOPT  .KQ.  2)  THEN 

0752 

,  DO  3220  I=1,P 

0753 

DO  3225  J=1,Q 

0754 

READ( 1 ,•)  PBAR(I, J) ,TSUBA(1, J) ,TAU(I, J) 

0755 

3225 

CONTINUE 

0756 

3220 

CONTINUE 

0757 

END  IF 

0758 

CLOSE(UNIT=1 ) 

0759 

* 

■RlNC=2.«MAXR/P 

0760 

0761 

CALL  BELL 

0762 

CALL  BELL 

0763 

CALL  BELL 

0764 

QRIDUP  =  .TRUE. 

0765 

GOTO  105 

0766 

0767 

END  IF 

0768 

0769 

0770 

c _ 

C  THE 

USER  CAN  WRITE  THE  CURRENT  GRID  TO  A  DATA  FILE  FOR  FUTURE 

0771 

0772 

0773 

C  USE 
c 

TO  AVOID  SPENi-ING  THE  TIKE  GENERATING  IT. 

0774 

3500 

CALL  SCLEAU(TRHNUH) 

0775 

WRITE(6,«) 

0776 

WRITE(6,#,‘  'Enter  the  name  you  wish  to  give  the  file.' 

0777 

READ( V'A20)')  FN IMS 

0778 

0PEN(UNIT=1 ,FILE=FNANE , STATUS* ' NEW ' ) 

0779 

WRITE(6,«) 

0780 

WRITE(6,«)  'Writing  grid  to  file  », FNAME 

0781 

0782 

0783 

0784 

0785 

0786 

0787 

0788 

0789 

0790 

0791 

0792 

0793 

0794 

0795 

0796 

0797 

0798 

0799 

0800 

0801 

0802 

0803 

0804 

0805 

0806 

0807 

0808 

0809 

0810 

0811 

0812 

0813 

0814 

0815 

0816 

0817 

0818 

0819 

0820 

0821 

0822 

0823 

0824 

0825 

0826 

0827 

0828 

0029 

0830 

0831 

0832 


WRITE (1,»)  IOPT 
WRITE(1,*)  MAXR 
IF  (IOPT  .EQ.  1)  THEN 
DO  3505  1=1, P 
DO  3510  J=1 ,Q 

WRITE(1,»)  PBAR(I| J) ,PREFL(X,J) ,T3GBA(X$J),TAU(I, J) 
3510  CONTINUE 

3505  CONTINUE 

END  IF 

IF  (IOPT  .EQ.  2)  THEN 
DO  3520  1*1, P 
DO  3530  Js1,Q 

WRITE ( 1  ,«)  PBAR(I,  J)  ,TSUBA(I,J)  jTAUd,J) 

3530  CONTINUE 

3520  CONTINUE 

END  IF 

CLOSE ( UNIT* 1) 

CALL  BELL 
CALL  BELL 
GOTO  105 


C  IN  SECTION  4000  THE  USER  PICKS  WHICH  QUANTITY  HE  WANTS  TO  PLOT. 

C  THE  GRID  USED  MAY  BE  ONE  GENERATED  DURING  THE  CURRENT  INTERACTIVE 
C  SESSION  OR  IT  MAY  BE  ONE  WHICH  SAVED  FROM  A  PREVIOUS  SESSION. 

C  THIS  FEATURE  ALLOWS  THE  USER  TO  RUN  THE  TIM-CONSUMING  PART 
C  OF  THE  PROGRAM  ON  A  NON-GRAPHICS  TERMINAL  AND  THEN  MOVE  TO  A 
C  GRAPHICS  TERMINAL  TO  QUICKLY  GENERATE  THE  CONTOUR  MAPS. 

C  THE  CONTOUR  LEVELS  TO  BE  PLOTTED  ARE  SPECIFIED 
C  BY  THE  USER  OR  DETERMINED  AUTOMATICALLY.  THE  CONTOURING  ALGOR- 
C  ITHM  THEN  SEARCHES  THROUGH  THE  GRID  TO  DETERMINE  THE  PATH  OF 
C  EACH  CONTOUR  LEVEL.  ONCE  THE  PATH  HAS  BEEN  DETERMINED,  IT  IS 
C  PLOTTED  USING  THE  DISSPLA  SUBROUTINES. 

4000  CALL  SCLEAR(TRMNUM) 

IF  (.NOT,  GRIDUP)  THEN 

WRITE(6,*)  'You  must  generate  or  read  a  grid  of  funotion' 
WRITE(6 ,*)  'values  before  plotting.  Hit  return  to  return' 
WRITE(6 ,*)  'to  main  menu.  Then  choose  menu  option  3* ' 
READ(V(F3.1)')  PAUS 
GOTO  105 
END  IF 

C - - - 

C  THE  LABELS  ARE  DEFINED. 

C - 


LABEL1= 'DISTANCE  FROM  ORIGIN  OF  PLANE,  aetera  • 

LABEL2= 'ANGLE  OF  ELEVATION  (deg) s' 

LABEL3=' HEIGHT  OF  MUZZLE  (meters)*' 

LABEL6='YOUR  MAXIMUM  CONTOUR  VALUE  IS  TOO  LARGE,  PLEASE 
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0833 
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0836 
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0838 
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0842 
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0844 
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0846 
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0849 

0850 

0851 

0852 

0853 

0854 

0855 

0856 

0857 

0858 

0859 

0860 

0861 

0862 

0863 

0864 

0865 

0866 

0867 

0868 

0869 

0870 

0871 

0872 

0873 

0874 

0875 

0876 

0877 

0878 

0879 

0880 

0881 

0882 

0883 

0884 


&  TRY  AGAIN.  IT  MUST  NOT  EXCEED 

LABEL7=»Y0UR  MAXIMUM  CONTOUR  VALUE  IS  TOO  SMALL,  PLEASE 
&  TRY  AGAIN.  IT  MUST  EXCEED 
LABEL8=» INPUT  NEW  VALUE  (ENTER  0.  FOR  DEFAULT).* 


NAME1=» CONTOURS  OF  PEAK  INCIDENT  OVERPRESSURE,  mbar* 
NAME2=  *  CONTOURS  OF  PEAK  REFLECTED  OVERPRESSURE,  mbar* 
NAME3=' CONTOURS  OF  BLAST  WAVE  TIME  OF  ARRIVAL,  mseo* 
KAME4=* CONTOURS  OF  POSITIVE  PHASE  DURATION,  msec* 


WRITE  (6,*) 
WRITE  (6,») 
WRITE  (6,*) 
WRITE  (6,») 
WRITE  (6,*) 
WRITE  (6 **> 
WRITE  (6,*) 
*  WRITE  (6,*) 
WRITE  (6,*) 
READ*, MENU 


•  PLOTTING  MENU* 

*  1.  Plot  contours.' 

'2.  Enlarge  or  reduce  plot  size.' 

'  3.  Deactivate/activate  polar  grid  on  plot.* 
'  4.  Return  to  proceeding  menu.' 

'Enter  the  corresponding  integer.' 


IF  (MENU  .EQ.  1)  GOTO  4100 
IF  (MENU  .EQ.  2)  GOTO  4200 
IF  (MENU  .EQ.  3)  GOTO  4300 
IF  (MENU  .EQ.  4)  GOTO  105 


C — - - - - - - - - - » - 

C  IF  THE  USER  WISHES  TO  PLOT  CONTOURS  ON  A  GRID  THAT  WAS  PREVIOUSLY 
C  SAVED  IN  A  DATA  FILE  HE  MAY  SPECIFY  THE  DATA  FILE  NAME. 

C - - - - - 

4100  CALL  SCLEAR(TRMNUM) 

C - - - - - - - 

C  THE  USER  CHOOSES  WHICH  BLAST  QUANTITY  HE  WANTS  TO  PLOT. 

C - - - 


Choose  which  quantity  you  wish  to  plot:' 


WRITE(6,«) 

WRITE(6,*) ' 

WRITE(6,») 

IF  (IOPT  .EQ.  1)  THEN 

WRITE(6,*)'  1.  Peak  incident  overpressure. * 

Peak  reflected  overpressure.* 

Blast  wave  time  of  arrival.' 

Blast  wave  positive  phase  duration. * 


WRITE(6,*) ' 
WRITE(6 ,*) ' 
WRITE(6,*) » 
WRITE(6,*) 
READ* , POPT 
ELSE 

WRITE(6,*)' 
WRITE(6,*)' 
WRITE(6 ,») ' 
WRITE(6,*) 


2. 

3. 

4. 


1.  Peak  incident  overpressure.' 

2.  Blast  wave  time  of  arrival.' 

3.  Blast  wave  positive  phase  duration.' 
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0891 
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0930 

0931 

0932 

0933 

0934 

0935 

0936 


READ*, POPT 

IF  (POPT  .EQ,  3)  POPT* 4 
IF  {POPT  .EQ.  2}  ?0PT=3 
END  IF 


C  THE  USER  CHOOSES  WHETHER  HE  WANTS  TO  SPECIFY  THE  VALUE  OP  EACH 
C  CONTOUR  OR  HAVE  THE  CONTOUR  LEVELS  AUTOMATICALLY  SCALED  STARTING 
C  WITH  A  COMPUTER  GENERATED  MINIMUM  VALUE. 

C - - — — - - - - - - - - - ...... _ _ _ 

CALL  SCLEAR(TRMNUM) 

WRITE(6f»)»You  may* 

WRITE(6,*) 

WRITE(6 ,*) *  1.  Enter  the  oontour  values  to  be  plotted.' 

WRITE(6 ,*) *  2.  Have  the  contour  values  automatically* 

WRITEC6,*)'  calculated  baaed  on  the  minimum  and  maximum* 
WRITE(6,*) •  values  of  the  function  in  the  plotting  domain.' 
WRITEC6,*) 

READ*,  MENU 
IF  (MENU.EQ.1)  THEN 


CALL  SCLEAR(TRMNUM) 

WRITE ( 6, *) 'Input  oontour  values  to  be  plotted  (mbar  or  mseo,' 
WRITE(6,*) 'max  of  40  values)' 

WRITE(6,*) 

WRITE (6,*) 'If  no  more  inputs,  enter  0.  ' 

4101  FORMAT ( '  », 'Contour  ',1X,I2) 

ICOUNT=0 
DO  4135  1=1,40 
WRITE(6 ,*) 

WRITE(6 ,4101 )  I 
READ  *,  H(I) 

IF  (POPT.EQ. 1 .AND.H(I) .GE.2.4)  THEN 
KLAPs 1 
ELSE 

KLAP=0 
END  IF 

H(I)  =  H(I)/1000. 

IF ( H( I ) .EQ.0.0)  THEN 
NLEVLSsICOUNT 
GO  TO  4136 
ELSE 

ICOUNT=IC0UNT+1 
END  IF 
4135  CONTINUE 


C - - - - - - - 

C  IF  THE  USER  INPUTS  THE  CONTOUR  VALUES  HIMSELF,  THEY  ARE 
C  SORTED  INTO  INCREASING  ORDER. 

C - - - 
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0937 

4136 

.  FINALsNLEVLS 

0938 

TESTsNLEVXS-1 

0939 

DO  4140  J=1,TEST 

0940 

LIMIT=FINAL-1  . 

0941 

DO  4145  1=1, LIMIT 

0942 

IF  (H(I).GT.H(I+1))  THEN 

0943 

TEMPsH(I) 

0944 

H(I)=H(I+1 } 

0945 

H(I+1 )=TEMP 

0946 

END  IF 

0947 

4145 

CONTINUE 

0948 

FINAL=FINAL-1 

0949 

4140 

CONTINUE 

0950 

0951 

ELSE 

0952 

0953 

C— — . 

0954 

C  IF  ' 

THE  CONTOUR  VALUES  ARE  TO  BE  COMPUTER  CHOSEN,  THE  MAX- 

0955 

C  IMUM  AND  MINIMUM  VALUES  IN  THE  CONTOUR  GRID  ARE  FOUND  AND 

0956 

C  THE 

CONTOUR  LEVELS  ARE  SCALED  ACCORDING  TO  THE  INCREMENT, 

0Q57 

C _ — . 

0958 

0959 

WRITE(6,«) 

0960 

WRITE(6,*)  'Enter  the  number  of  contours  you  want  plotted 

0961 

WRITE(6,#)  '(maximum  of  40).' 

0962 

READ*, NLEVLS 

0963 

WRITE(6,») 

0964 

0965 

4250 

IF  .(P0PT.EQ.1)  THEN 

0966 

CALL  FNDMAX ( P8AR , PMAX , PO  PT ) 

0967 

CALL  FNDMIN ( PBAR , PMIN , POPT ) 

0968 

PMIN  =  TWOSIG( 1 , 1*PMIN) 

0969 

PMAX  =  ,9*PMAX 

0970 

CINC= ( PMAX-PMIN ) / ( NLEVLS- 1 ) 

0971 

CINC=ONESIG(CINC) 

0972 

DO  4255  1=1, NLEVLS 

0973 

H(I)=PMIN  +  ( I— 1 ) *CINC 

0974 

IF  (H(I)  ,GT.  2.4)  THEN 

0975 

FLAP  =  1 

0976 

ELSE 

0977 

KLAP  =  0 

0978 

END  IF 

0979 

4255 

CONTINUE 

0980 

END  IF 

0981 

IF  (P0PT.EQ.2)  THEN 

0982 

CALL  FNDMAX ( PREFL , PRMAX , POPT ) 

0983 

CALL  FNDMIN ( PREFL, PRMIN, POPT) 

0984 

PRMIN  =  TW0SIG( 1 . 1*PRMIN) 

0985 

PRMAX  =  0 . 9#PRMAX 

0986 

CINC=( PRMAX-PRMIN ) / ( NLEVLS- 1 ) 

098? 

CINCsONESIG ( CINC ) 

0988 

DO  4260  1=1, NLEVLS 
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0989 

0990 

0991 

0992 

0993 

0994 

0995 

0996 

0997 

0998 

0999 

1000 

1001 

1002 

1003 

1004 

1005 

1006 

1007 

1008 

1009 

1010 
1011 
1012 

1013 

1014 

1015 

1016 

1017 

1018 

1019 

1020 
1021 
1022 

1023 

1024 

1025 

1026 

1027 

1028 

1029 

1030 

1031 

1032 

1033 

1034 

1035 
103' 

1037 

1038 

1039 

1040 


H(I)=PRMIN  +  ((I-1)*CINC) 

4260  CONTINUE 
END  IF 

IF  (P0PT.EQ.3)  THEN 

CALL  FNDMAX ( TSUB A , TSAMAX , FQPT ) 

CALL  FNDMIN ( TSUBA , TSAMIN , POPT ) 

TSAMIN  =  TSAMIN*  1.1 
TSAMAX  »  TW0SIG( TSAMAX*. 9) 

CINCs ( TSAMAX- TSAMIN } / ( NLEVLS- 1 ) 

CINCaONESIQ ( CINC ) 

DO  4265  1-1 ,NLEVLS 

H(I) sTSAMAX-( (1-1 ) *CINC) 

4265  CONTINUE 
END  IF 

IF  (P0PT.BQ.4)  THEN 

CALL  FNDMAX (TAU,TAUMAX, POPT) 

CALL  FNDMIN ( TAU , TAUMIN , POPT ) 

TAUMIN=TAUMIN* 1 , 1 
TAUKAXsTWOSIu ( TAUMAX* . 9) 

CINCs  ( TAUMAX-TAUMIN )  /  ( MLEVLS- 1 ) 

CINCsONESIG(CINC) 

DO  4267  1=1 iNLEVLS 

H ( I ) sTAUMAX- ( ( I- 1 ) *CINC ) 

426?  CONTINUE 
END  IF 

END  IF 

C  FOUlw(I)=0  MEANS  THAT  NO  CONTOURS  WERE  LOCATED  AT  THE  H(I) 

C  LEVEL.  FOUND(I)o1  MEANS  THAT  CONTOURS  HERE  LOCATED. 

DO  4173  1=1 ,NLEVLS 
FOUND(I)=0 
4173  CONTINUE 

C - - - — - — - - 

C  THE  USER'S  TERMINAL  IS  READIED  FOR  THE  PLOTS  .THE  CONTOURS 
C  WILL  BE  DRAWN  USING  DISSPLA  VERSION  9.0  COMMANDS.  BECAUSE 
C  OF  THE  SYMMETRY  OF  THE  PROBLEM,  THETA  NEED  ONLY  RANGE  FROM 
C  0  TO  180  DEGREES. 

IF  (TRMNUM.EQ.O  .OR.  TRMNUN  .EO.  4)  THEN 
CALL  RETRO  (BAUD) 

ELSE  IF  (TRMNUM.EQ. 1 )  THEN 
CALL  TK4010  (BAUD) 

ELSE  IF  (TRMNUM.EQ. 2)  THEN 
CALL  TK4014  (BAUD,1) 

ELSE 

CALL  HP2623  (BAUD) 
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END  IP 
XSTEPsMAXR/5 


1041 

1042 

1043 

1044 

1045 

1046 
104? 

1048 

1049 

1050 

1051 

1052 

1053 

1054 

1055 

1056 

1057 

1058 

1059 

1060 
1061 
1062 

1063 

1064 

1065 

1066 

1067 

1068 

1069 

1070 

1071 

1072 

1073 

1074 

1075 

1076 

1077 

1078 

1079 

1080 
1081 
1082 

1083 

1084 

1085 
1036 

1087 

1088 

1089 

1090 

1091 

1092 


C - - - - - - 

C  THE  DISSPLA  COMMANDS  READY  THE  TERMINAL  FOR  PLOTTING 
C  BY  FIRST  DEFINING  THE  SUBPLOT  AREA  THEN  CALIBRATING  IT 
C  FOR  THE  CONTOUR  LINES. 

C - - - - - „ - 


CALL  BGNPL(O) 

CALL  NOCHEK 
CALL  NOBRDR 
CALL  SIMPLY 
CALL  PHYSOR(. 1,1.0) 

CALL  BLOWUP (XYF AC) 

CALL  TITLE ( '  $',-100,'  « ,0, •  ',0,10,0,5.0) 
CALL  FRAME 

CALL  GRAF(0.0,XSTEP,2.»MAXR,0.OfXSTEP,MAXR) 


C - - - - - - - — — - — — _ 

C  TH7  CONTOUR  SEARCH  SUBROUTINE  IS  CALLED  UITH  THE  APPRO- 
C  PRIATR  QUANTITY  TO  BE  PLOTTED. 


IF  (POPT.EQ. 1 )  THEN 

CALL  CONTOR(PBAR, POUND) 
END  IF 

IF  (POPT.EQ. 2)  THEN 

CALL  CONTGR ( PREFL , FOUND ) 
END  IF 

IF  (POPT.EQ. 3)  THEN 

CALL  CONTORfTSUBA, FOUND) 
END  IF  . 

IF  (POPT.EQ. 4)  THEN 

CxVLL  CONTOR(TAU, FOUND) 
END  IF 


C - - - - - — ... - 

C  THE  APPROPRIATE  HEADING  IS  PLACED  ON  THE  PLOT. 


IF  (POPT.EQ.1)  THEN 
HCHAR=44 
NAHF=NAK£ 1 
END  IF 

IF  (POPT.EQ. 2)  THEN 
;<CHAS-45 
NAHE=NakS2 
END  IP 

IF  (POPT.EQ. 3)  TyEN 
NCHARs44 
NAHE=NAME3 


1093 

1094 

1095 

1096 

1097 

1098 

1099 

1100 
1101 
1102 

1103 

1104 
1  i  05 
1106 

1107 

1108 

1109 

1110 
1111 
1112 

1113 

1114 

1115 

1116 
1  M? 
1118 

1119 

1120 
1121 
1122 

1123 

1124 

1125 

1126 
112? 
1128 

1129 

1130 

1131 

1132 

1133 

1134 

1135 

1136 

1137 

1138 

1139 

1140 

1141 

1142 

1143 

1144 


END  IF 

IF  (POPT.EQ.4)  THEN 
NCHAR=42 
NAME=NAME4 
END  IF 

CALL  HEIGHT (0.22) 

CALL  MESSAG(*REF(NAME),NCHAR,.75,7.2) 
CALL  HEIGHT (0.19) 


C - — - - 

C  THE  REMAINING  PLOT  SPECIFICATIONS  ARE  WRITTEN  ON  THE  PLOT. 
C - - - - - - 


CALL  MESSAG(f REF(GON) ,32,4.5,6.8) 

CALL  HESSAG(*REF(LABEL2)l25f4.5,6.5) 

CALL  RE ALNO ( f REF ( ELEANG ), 2, 8. 6,6. 5) 

CALL  HESSAG(*REF(LABEL3),26,4.5,6.2) 

CALL  REALNQ ( *REF ( ZHEIGH ) , 2 , 8 . 6 , 6 . 2 ) 

CALL  MESSAG( * ASUBM(n/s )  =  • ,  1 1 , 0 . 0 , 6 .8) 

CALL  REALNO( *  REF ( ASUBM) , 1 , 2 . 0 , 6 . 8 ) 

CALL  HESSAG( 'PRATIO=*,7,O.Q,6.5) 

CALL  REALN0(*REF(PRATI0),1 ,2. 0,6. 5) 

CALL  MESSAG( 'VSUBP(o/3)=» ,1 1 ,0.0,6. 2) 

CALL  REALNO($REF(VSUBP) , 1 ,2.0,6 .2) 

CALL  MESSAG( 'GAMMA=* ,6,0. 0,5. 9) 

CALL  REALNO(*REF(GAHHP),2,2.0,5.9) 

CALL  MESS AG ( 'MU =  ',6,0. 0,5. 6) 

CALL  REALNO(f REF(MU) ,2, 2. 0,5. 6) 

IF  (INTBAL)  THEN 

CALL  MBSSAG( 'PROPELLANT  MASS(ks)e ' ,20,4.5,5.9) 

CALL  REALNO{ $ R8F ( PROPM) , 5 , 7 . 7 , 5 . 9 ) 

CALL  KESSAG{ 'PROJECTILE  HASS(kg)s’ ,20,4.5,5.6) 

CALL  REALNO(*REF(PRQJM),4,7.7,5.6) 

CALL  MESSAG( 'SPECIFIC  FORCE(  )a» ,24,0.0,5.2) 

CALL  MIXALF( 'INSTRUCTION' ) 

CALL  HESSAG(  ,ia(EH.6)2(EXHX)/8(EH.6)2(MHX)$' ,  100, 

A  2.00,5.2) 

CALL  RBSET< 'MXXALF* ) 

CALL  REALNO( *REF ( SPSCF ), 104, 3. 20, 5. 2) 

END  IF 

CALL  MESSAGE  REF  ( LABEL  1 ) ,38, 2. 25, -.75) 

CALL  ENOGR(O) 

CALL  HEIGHT( .25) 

IF  (NOQRID.EQ.O)  THEN 
CALL  Y NON OH 

CALL  TITLE ( •  $',-100,'  * ,  1 , *  *,1,10.0,5.0) 

CALL  POLAR ( DTOR , XSTEP , 5 . 0 , 0 . 0 ) 

GO  TO  4122 
END  IF 

CALL  TITLEC  $',-100,'  ',1,*  *,0,10.0,5.0) 
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.  i  V"  A,  s  »  «  *■ \  A' ■  ‘•T-_  A,  A-  '  » 


*  O'  fJl  C  V**  O. ’ 


1145  C- - — - — — - - - - 

1146  C  THE  ACTUAL  CONTOUR  LINES  ARE  DRAWN  IN  THE  LINEAR  MODE. 

1147  C  NOW,  THE  POLAR  MODE  IS  CALLED  IN  ORDER  TO  SUPERIMPOSE 

1148  C  THE  DISTANCE  SCALE  AND  ADD  THE  GRID  IF  ACTIVATED. 

1 1 49  C - — — - - - - - — - - - — 

1150 

1151  CALL  POLAR ( DTOR , XSTEP , 5 . 0 , 0 . 0 ) 

1152  RATI0(1)c1. 

1153  RATI0(2)=1 . 

1154  CALL  MRSCOD ( » 3 » 2 , RATIO ) 

1155  CALL  GRID (-30, 1 ) 

1156  4122  CALL  ENDPL(O) 

1157 

1158 

1159 

1160 
1161 
1162 

1163 

1164 

1165 

1166  C 

1167  C  TO ERASE  THE  PLOT  AND  RETURN  TO  PROGRAM  EXEUTION,  THE  USER 

1168  C  MUST  PRESS  'RETURN*. 

1169  C— , - - - - - - - - - - - 

1170 

1171  WRITE(6,*} 

1172  CALL  SCLEAR(TRMNUM) 

1173  IF  (POPT.SQ.1 .AND.KLAP.EQ.1 )  THEN 

1174  WRITE(6 ,*) 

1175  WR1TE(6 ,*) ’The  scaling  relations  used  in  this  program  were' 

1176  WRITE (6,*) 'derived  using  pressure  data  of  no  more  than' 

1177  WR1TE(6,*) ’2.4  atm.  However,  this  graph  exceeds  this  pressure' 

1178  WRXTE(6 ,*) 'and  these  higher  pressure  predictions  may  be' 

1179  WRITE (6,*) 'considered  less  aoourate.* 

1180  WRITE(6 ,*) 

1181  WRITE(6,*) 'Press  return  to  oontinue. ' 

1182  READC5,*)  DUMMY 

1183  END.  IF 

1184  IF  (P0PT.EQ.3.AND.KLAP.EQ.1)  THEN 

1185  WRITE(6 ,*) 

1186  WRITE(6,*) 'Some  of  the  contours  exceed  the  data  range  used* 

1187  WR1TE(6 ,*) 'in  deriving  the  soaling  relations.  Therefore,' 

1188  WRITE(6 ,*) 'The  lower  time  values  may  be  less  aoourate* 

1189  WRITE(6 ,*) 

1190  WRITE(6,*) 'Press  return  to  continue.* 

1191  .  READ(5, ' (A1 ) »)  DUMMY 

1192  END  IF 

1193 

1194  GOTO  4000 

1195 

1196  4200  CALL  SCLEAR(TRMNUM) 


IF  (TRMNOM.EQ.O)  CALL  XRETRO  (480) 
IF  (TRMNUM.EQ.3)  PRINT*, ESC// »*dA« 
IF  (TBMNUM  .EQ.  4)  THEN 
CALL  XRE!TR0(480) 

PRINT*, ESC// *2' 

END  IF 
CALL  DONEPL 


1197 

1198  C — - 

1199  C  THE  PLOT  SIZE  CAN  BE  ALTRERED. 

1 200  C - - - - - 

1201 

1202  WRITE(6,*) 

1203  WRITE(6, •) 

1204  WRITE(6,*) 'Input  faotor  for  ohaaging  graph  size.' 

1205  WRITE (6,*) 'Enter  1.0  to  obtain  the  default  size.' 

1206  WRITE(6,«) 

1207  READ  *,XYFAC 

1208  GOTO  4000 

1209 

1210  4300  CALL  SCLEAR(TRMNUH) 

1211 

1212  C - - - „ - - * _ 

1213  C  THE  USER  CAN  ELIMINATE  OR  REACTIVATE  THE  POLAR  GRID  THAT  IS 

1214  C  SUPERIMPOSED  ON  THE  FLOT. 

1215  C - - - - - - - 

1216 

1217  WRITE(6,*) 

1218  WRITE(6,a) 

1219  WRITE(6,*)'A  polar  grid  is  by  default  superimposed  on  the' 

1220  WRITE(6,*) 'oontour  maps.  To  erase  this  grid,  enter  OFF,' 

1221  WRITE (6,*) 

1222  WRITE(6,*) 'If  the  grid  is  already  deactivated,  enter  ON,’ 

1223  WRITE(6, •) 'Hit  return  for  no  ohange.* 

1224  WRITE(6,*>) 

1225  READ(5,'(A3)')  ERASE 

1226  IF  ( ERASE. EQ. 'OFF')  NOGRIDsO 

1227  IF  (ERASE. EQ. 'OH  »)  NOGRID a 1 

1228  GO  TO  4000 

1229 

1 230  C - „ - * - - 

1231  C  SECTION  5000  WRITES  THE  INPUT  INFORMATION  TO  A  PILE  FOR 

1232  C  FUTURE  USE. 

1234 

1235  5000  CALL  SCLEAR(TRMNUM) 

1236  WRITE(6,»)  'Enter  the  name  you  wish  to  give  the  file.' 

1237  READ(VU15)'}  WFILE 

1238 

1239  OPEN  (UNIT=1  ,FILEsWFXLE, STATUS* 'NEW* ) 

1240  WRITE (1,'(A32>»)  GUN 

1241  WRITE(1,*)  D 

1242  WRITE( 1 , * ( A 1 ) » )  ANSI 

1243  IF  (ANSI  .UQ.  'Y* )  THEN 

1344  WRITE( 1 ,*)  PROPM 

1245  WRITE( 1 ,•)  TRAVEL 

1246  WRITE ( 1 ,•)  CHAMV 

124?  WRITE(1,*)  PROJ 

1248  WRITE( 1 , •)  PROPB 
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1249 

WRITE (1 , » (A1 ) ® )  WUNITS 

1250 

WRITEO  ,*)  SPECF 

1251 

WRITEO,*)  TEMPK 

1252 

ELSE 

1253 

WRITEO,*)  ASUBM 

1254 

WRITEO,*)  PRATIO 

1255 

END  IF 

1256 

WRITEO,*)  AM3SS 

1257 

WRITEO,*)  VSUBP 

1258 

WRITEO,*)  GAMMP 

1259 

WRITEO,*)  MU 

1260 

WRITEO,*)  ELEANG 

1261 

WRITEO ,»)  ZHEIGH 

1262 

CLOSE(UNIT=1 ) 

1263 

1264 

GOTO  105 

1265 

1266 

6000  CONTINUE 

1267 

1268 

STOP 

1269 

END 

0001 

0002  C - ....................... 

0003  e  SUBROUTINE  CONTOR  SEARCHES  FOR  THE  fc&GINNXNO  OF  CONTOURS  BY  COH- 
0004  C  PARING  THE  VALUE  OF  THE  FUNCTION  AT  ADJACENT  GRID  SQUARES.  IT  IS 
0005  C  CALLED  BY  THE  MAIN  PROGRAM  ONCE  THE  FUNCTIONAL  VALUES  AT  ALL  THE 

0006  C  GRIDPOINTS  HAVE  BEEN  DETERMINED. 

0007  C - - - - - 

0008 

0009  SUBROUTINE  CONTOR(GRIDPT, FOUND) 

0010  PARAKETER(NL=40,Qs50,Pb100) 

0011  INTEGER  XM(P)  ,YM(Q)  ,UNUSED(P,Q)  ,FWND(NL)  ,OPEN 

0012  INTEGER  P0PT,PP,QQ 

0013  REAL  H(NL),QRIDPT(FfQ),PBAR(P,Q),Ht<R 

0014  COMMON/ GRID/  NLEVLS , H , & ,  OPEN , XM , YH ,  .UtXQN 

0015  COMMON/ INC/ RING ,MAXR 

0016  COMMON/ PRESS/ PBAR , POPT  , KL AP 

0017  COMMON  WQRK( 12000) ,XPL0T(250) ,Y PLOT (250) , 

0018  $  LSAV , KSAV , NOSTOR ( NL ) ,HSAV 

0019 

0020  HFACT  a  1000.0 

0021 

0022  PP  a  P 

0023  QQ  =  Q 

0024 

0025  &0  VI  r  1*1  ,.P 

0026  490  mm 

0027  DO  »20  I»1,Q 
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. 


0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0070 

0079 


420  YM(I)=I 
JMsQ-1 
IM=P-1 

C - - - - - - - - - 

C  ALL  CONTOURS  AT  EACH  LEVEL  ARE  DRAWN  BEFORE  MOVING  ON  TO  THE  NEXT 
C  LEVEL. 

C- - - - - - 

C - - - - - - - - - 

C  THE  ’ XPLOT '  AND  T YPLOT1  VALUES  WHICH  HAVE  BEEN  STORED  IN 
C  THEIR  RESPECTIVE  ARRAYS  DURING  SUBROUTINE  PLOT  ARE  FED 
C  INTO  DISSPLA’S  CONTOUR  SETUP  OPTION  *CONCRV>,  ALONG  WITH 
C  THE  APPROPRIATE  PRESSURE  LEVEL. 

C - — - - - - 

CALL  BCOMON( 12000) 

CALL  CONBGN 
MSAV  =  1 

DO  430  K=1,NLEVLS 

KSAV=K 

LSAV=0 


C - - - - - 

C  THE  ARRAY  UNUSED  IS  INITIALIZED  FOR  THE  CONTOUR  LEVEL.  UNUSED ( I, J)= 
C  1  MEANS  THAT  THE  I, J  GRID POINT  WILL  BE  USED  IN  DRAWING  THE  CONTOURS 
C  AT  THIS  LEVEL.  UNUSED(I, J)=0  MEANS  THAT  THE  I,J  GRID POINT  IS  NOT 
C  SIGNIFICANT.  .  AT  EACH  CONTOUR  LEVEL,  THE  BOUNDARY  OF  THE  GRID  IS 
C  SCANNED  FOR  A  POINT  WHERE  A  CONTOUR  CROSSES  INTO  THE  GRID  WITH  THE 
C  HIGHER  FUNCTIONAL  VALUES  ON  THE  RIGHT.  THIS  CONDITION  PREVENTS  THE 
C  PROGRAM  FROM  RELOCATING  CONTOURS  IT  HAS  ALREADY  DRAWN.  THE  SCAN  OF 
C  THE  GRID  BOUNDARY  BEGINS  IN  THE  LOWER  LEFT  CORNER  AND  PROCEEDS  COUN- 
C  JTERCLOCKWISE , 

C 

C  IF  AN  INTERSECTION  WITH  THE  GRID  BOUNDARY  IS  FOUND,  SUBROUTINE  FOL- 
C  LOW  IS  CALLED  TO  FOLLOW  THIS  CONTOUR  THROUGH  THE  GRID  UNTIL  IT 
C  EXITS  FROM  THE  GRID. 

C - 

DO  440  J-2, JM 
DO  445  1=2, IM 
UNUSED ( I, J)=0 

IF  (GRIDPT(I-1 ,  J) .LT.H(K) . AND, GRID PTC  I, J) .GE.H(K) ) 

A  UNUSED(I,J)=1 

445  CONTINUE 

440  CONTINUE 
OPEN=1 


C - - - 

C  FIRST,  THE  BOTTOM  EDGE  OF  THE  GRID  IS  SCANNED. 
C - 


DO  450  1=2, P 
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VO 


1 1 


0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 


IF  <GRIDPT(I-1,1).LT.H(K).AND.GRIDPT(I,1).GE.H(K))  THEN 
F0UND(K)=1 

CALL  FOLLOW ( GRIDPT , I , 1 , -1 , 0 , UNUSED) 

CALL  CONCRV(XPLOT(1) ,YPL0T(1) ,LSAV,HFACT*H(K)) 

LSAV  =  0 
END  IF 
450  CONTINUE 

C - 

C  THE  RIGHT  EDGE  IS  SCANNED. 

C - 

DO  460  J=2,Q 

IF  (GRIDPT (P,J-1 ) .LT.H(K) .AND.GRIDPT(P, J) .GE.H(K))  THEN 
F0UND(K)=1 

CALL  FOLLOW ( GRIDPT , PP , J , 0 , - 1 , UNUSED ) 

CALL  CONCRV(XPLOT( 1 ) , YPL0T( 1 ) , LSAV , HFACT*H ( K ) ) 

LSAV  =  0 
END  IF 

460  CONTINUE 

C - - - * - 

C  THE  TOP  EDGE  IS  SCANNED. 

C - - - - 

DO  470  L=1,IM 
IsP-L 

IF  (GRIDPT(I+1 ,Q) .LT.H(K) .AND. GRID PT(I,Q) .GE.H(K))  THEN 
FOUND(K)=1 

CALL  FOLLOW ( GRIDPT , I , QQ , 1 , 0 , UNUSED ) 

CALL  CONCRV(XPLOT( 1 ) , YPLOT( 1 ) ,LSAV,HFACT*H(K) ) 

LSAV  =  0 
END  IF 
470  CONTINUE 

C - 

C  THE  LEFT  EDGE  IS  SCANNED. 

C - 

DO  480  L=1,JN 
J=Q-L 

IF  ( GRIDPT ( 1 ,J*1 ) .LT.K(K) . AND.GRIDPT( 1,J).GE.H(K))  THEN 
FOUND(K)s1 

CALL  FOLLOW ( GRIDPT, 1,J, 0,1, UNUSED) 

CALL  CONCRV(XPLOT( 1 ) ,YPL0T( 1 ) ,LSAV,HFACT*H(K) ) 

LSAV  =  0 
END  IF 

480  CONTINUE 

- - 

'  C  ONCE  ALL  THE  OPEN  CONTOURS  AT  LEVEL  H  ARE  FOUND  AND  FOLLOWED  TO 
C  COMPLETION,  THE  ARRAY  UNUSED  IS  SEARCHED  FOR  ANY  CLOSED  CONTOURS 
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0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 


C  OF  HEIGHT  H.  IF  A  CLOSED  CONTOUR  IS  FOUND,  SUBROUTINE  FOLLOW  IS 
C  CALLED  TO  FOLLOW  IT  UNTIL  IT  RETURNS  TO  THE  POINT  WHERE  IT  WAS 
C  FOUND. 

C - - - - - 

MSAV  =  2 
0PEN=0 

DO  490  L=2, JM 
J=Q-L+1 
DO  500  M=2,IM 
I=P-M+1 

IF  (UNUSED(I.J).EQ.I)  THEN 
FOUND (K) =1 

CALL  FOLLOW ( GRIDPT , I , J , - 1 , 0 , UNUSED ) 

CALL  CONCRV(XPLOT( 1) ,YPL0T(1 ) ,LSAV,HFACT*B(K)) 
LSAV  =  0 
END  IF 

500  CONTINUE 

490  CONTINUE 
430  CONTINUE 

CALL  CONEND 

C  NEXT,  THE  DISSPLA  PARAMETERS  ARE  SET  TO  DRAW  SOLID  CONTOUR 
C  LINES  AND  LABEL  THEM  WHERE  POSSIBLE.  THE  *CONTUR'  COMMAND 
C  ACTUALLY  DRAWS  THE  LINES  USING  THE  DATA  STORED  IN  'COHCRY*. 


CALL  HEIGHT (.16) 

CALL  CONLIN ( 0,' SOLID',* LABELS »,1, 4) 
CALL  CONANG(90.) 

CALL  CONMIN (4.5) 

CALL  CONTUR(  1 ,  'LABELS* ,  'DRAW* ) 

RETURN 

END 


0001 

0002  C - - - — - — - — - 

0003  C  SUBROUTINE  FOLLOW  FOLLOWS  A  CONTOUR  THROUGH  THE  GRID  ONCE  IT  HAS 

0004  C  BEEN  FOUND  BY  CHECKING  THE  SIDES  OF  THE  GRID  SQUARES  TO  DETERMINE 

0005  C  WHICH  SIDE  THE  CONTOUR  PASSES  THROUGH.  IT  FINDS  THE  POINT  OF 

0006  C  INTERSECTION  USING  LINEAR  INTERPOLATION  BBTHBKN  THE  ADJACENT  GRID- 
0007  C  POINTS.  CONTROL  DOES  NOT  RETURN  TO  SUBROUTINE  CONTOR  UNTIL 

0008  C  THE  CONTOUR  EITHER  EXITS  FROM  THE  GRID  OR  RETURNS  TO  THE  POINT 

0009  C  IN  THE  GRID  WHERE  IT  WAS  FIRST  FOUND.  WHEN  SUBROUTINE  FOLLOW 

0010  C  DETERMINES  HOW  THE  CONTOUR  PASSES  THROUGH  A  GRID  SQUARE,  SUBROU- 

001 1  C  TINE  PLOT  IS  CALLED  TO  STORE  THE  X  AND  Y  VALUES  LATER  USED  FOR 

0012  C  THE  ACTUAL  PLOTTING. 


0015  SUBROUTINE  FOLLOW ( GRID PT, I IG, JJG.IIA, JJA.UNUSED) 

0016  PARAMETER ( NL=40 , Q=  50 , P= 1 00 ) 

0017  INTEGER  IG, JG,IA, JA,OPEN,K, LAST, FIRST, TE,XM(P) ,YM(Q) , 

0018  &  NLEVLS, UNUSED (P,Q) 

0019  REAL  GRIDPT{P,Q),H(NL),Z,ZA,ZB,ZC,X,Y, 

0020  &  RINC,T,MAXR 

0021  COMMON/ GRID/NLEVLS ,H,K, OPEN , XM , YM , ISIGN 

0022  COMMON/ INC/RINC,MAXR 

0023  COMMON  WORK( 12000), XPL0T(250),YPLOT(250), 

0024  A  LSAY , KSAV , NOSTOR ( NL ) , MSAV 

0025 

0026  C - - - - 

0027  C  THE  PARAMETERS  IG, JG,IA,  AND  JA  ARE  USED  IN  VARIOUS  COMBINATIONS 

0028  C  AS  THE  SUBSCRIPTS  FOR  THE  ARRAYS  XM  AND  YM.  THE  VALUES  OF  THESE 

0029  C  SUBSCRIPTS  DETERMINE  WHICH  GRID  POINTS  ARE  CONSIDERED  IN  ANY  PARTI- 

0030  C  CULAR  LINEAR  INTERPOLATION. 

0031  C - 

0032 

0033  IG=IIG 

0034  JGsJJG 

0035  IAeIIA 

0036  JA=JJA 

0037 

0038  C - 

0039  C  FIRST= 1  INDICATES  THAT  A  NEW  CONTOUR  HAS  JUST  BEEN  FOUND;  FIRSTsO 

0040  C  MEANS  THAT  THE  PLOTTING  OF  THE  CONTOUR  IS  ALREADY  IN  PROGRESS. 

0041  C  LAST= 1  INDICATES  THAT  THE  END  OF  A  CONTOUR  HAS  BEEN  REACHED;  LASTsO 

0042  C  MEANS  THAT  THE  CONTOUR  IS  NOT  YET  FINISHED. 

0043  C - - - - 

0044 

0045  FIRSTs 1 

0046  LASTsO 

0047 

0048  C— - - - - - 

0049  C  THE  LOCATION  OP  THE  POINT  T  WHERE  THE  CONTOUR  PASSES  BETWEEN  THE 

0050  C  POINTS  (I,J)  AND  (I+IA, J+JA)  IS  CALCULATED  USING  LINEAR  INTERPOLATION. 

0051  C— 

0052 
0053 
0054 

0055  605 

0056 
0057 
0058 
0059 

0060  C - 

0061  C  TESTS  ARE  NOW  PERFORMED  TO  DETERMINE  IP  T  IS  THE  LAST  POINT  ON  THE 
0062  C  CONTOUR. 

0063  C - - - 

0064 

0065  IF  ( OPEN , EQ . 1 )  GOTO  610 

0066  IF  (IA.EQ.-1 . AND. UNUSED (IG,JG) .EQ.O)  LASTs 1 
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ZsQRIDPT(IG , JG) 
ZAsGRIDPTCIG+IA, JG+JA) 

T=0. 

IF(Z.NE.ZA)  T=(Z-H(K))/(Z-ZA) 
XaXH(IG)-T*(XH(IG)-XM(IG+IA)} 
YsYN( JG)-T*(YH( JG)-YH(JG+JA) ) 


0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 

0117 


GOTO  650 

IP  (FIRST.EQ.1)  GOTO  669 
IP  (JA.EQ.O)  GOTO  620 
GOTO  630 

IF  (JG.EQ.1 .OR.JG.EQ.Q)  LASTs 1 
IP  (JA.NE.O)  GOTO  640 
GOTO  650 

IF  (IG.EQ.1 .OR.XG.EQ.F)  LASTsI 
IF  (LAST.EQ.1)  GOTO  660 
IP  (IA.EQ.-1 )  UNUSED ( IG , JG ) sO 


C  THE  COORDINATES  OF  T  ARE  OUTPUT  TO  SUBROUTINE  PLOT. 
C - - - — — - 


LSAV=LSAV+1 

CALL  PLOT(X,Y, FIRST) 


C  COMPARISONS  ARE  NOV  MADE  TO  DETERMINE  WHICH  OF  THE  CELL  SIDES 
C  THE  CONTOUR  CROSSES  NEXT.  THE  VALUES  OF  Z,  ZA,  10,  JO,  Ik, 

C  AND  JA  ARE  ADJUSTED  BEFORE  CONTINUING  TO  FIND  A  NEW  POINT,  T. 

IF  (LAST.EQ.1)  THEN 
NOS'i'OR  ( KSAV )  sLSAV 
RETURN 
END  IF 

ZBsGRIDPT ( IG+ JA , JG-XA) 

IF  (ZB.GE.H(K))  GOTO  67G 

ZAsZB 

TEsIA 

IA=JA 

JA=-TE 

GOTO  690 

670  ZC*GRIDPT(IG+IA*JA,JG-IA+JA) 

IF  (ZC.GE.H(K) )  GOTO  680 
Z-ZB 
ZAaZC 
IGnIG+JA 
JQaJG-IA 
GOTO  690 
680  Z=ZC 

IGaIG+IA+JA 
JGaJG-IA+JA 
TEsJA 
JAsIA 
IAs-TE 
690  FIRSTsO 
GOTO  605 
END 
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0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 


C  FUNCTIONS  FOR  WHICH  CONTOURS  ARE  PLOTTED. 

C - - - - - - 

SUBROUTINE  FNDMIN ( VALS2 , MINVAL , POPT ) 

PARAMETER ( NL=40 , Q=50 , P= 1 00 ) 

INTEGER  POPT 

REAL  VALS2(P,Q), MINVAL 

MINVALr 100000. 

DO  810  1=1 ,P 
DO  8 20  J=1 |Q 

IF(POPT.EQ.3.AND.I.GE.45.AND.I.LE.55,AND.J.LE.5) 
&  GO  TO  820 

IF(POPT.EQ.4.AND. I.GE.45.AND.I.LE.55.AND, J.LE.5) 
&  GO  TO  820 

IF  (VALS2(I,J).LT. MINVAL)  MINVALaVALS2(I, J) 

820  CONTINUE 
81 0  CONTINUE 
RETURN 
END 


0001 

0002  C - — - - - - 

0003  C  SUBROUTINE  REFCAL  IS  CALLED  BY  THE  MAIN  PROGRAM  TO  DETERMINE 

0004  C  THE  REFLECTED  OVERPRESSURE,  PREREF  AT  EACH  GRIDPOINT  IN  THE  CONTOUR 

0005  C  GRID.  FIRST  A  CHECK  IS  MADE  TO  SEE  IF  REGULAR  REFLECTION  IS 

0006  C  POSSIBLE.  IF  IT  IS  POSSIBLE,  THE  REGULAR  REFLECTION  CALCULATION 

0007  C  PROCEEDS.  IF  REGULAR  REFLECTION  DOES  NOT  OCCUR  SUBROUTINE 

0008  C  NONREG  IS  CALLED  TO  DETERMINE  THE  REFLECTED  OVERPRESSURE. 

0009  C - - - - - - - 

0010 

0011  SUBROUTINE  REFCAL ( PRESIN, PREREF, ALPHA 1 ) 

0012  REAL  MACHN1 ,MACHN2,LOVBND, PSIANG( 1 1 ) , ANGLIN ( 1 1 ) ,5L0PES( 11) 

0013  REAL  PSISLP(II) 

0014  COMMON/ REFL/ GAMMA , TOL , UPBND , LOWBND , XI 

0015  COMMON/ TRIG/ RTOD,DTOR, PI 

0016  COMMON/ LAG RAN / PSIANG , PSISLP , ANGLIN, SLOPES 

0017  CALL  PRELIM(ALPHA1 , PRESIN, DELTA, MACHN1 ,MACHN2,DELMAX) 

0018  IF  ( DELTA. GE.0.90«DELMAX)  THEN 

0019  CALL  NONREG (DELMAX, PRESIN, ALPHA1 ,MACHN1 .MACHN2, PREREF) 

0020  RETURN 

0021  END  IF 

0022  AL  PMAX  sF ALMAX ( MACHN2 , GAMMA ) 

0023  CALL  ITRATE(0,ALPHA2, DELTA, LOWBND, ALPMAX,HACHN2, PRESIN) 

0024  QUOTNTs ( 2 . •GAMMA* ( MACHN2t#2* ( SIN ( ALPHAS ) ) **2-1 . ) ) / 

0025  A  (GAMMA+I.HI . 

0026  PREREF=QU0TNT*XI-1 . 

0027  RETURN 

0028  END 


li 

,1 


0001 

0002  C. 

0003  C  SUBROUTINE  PRELIM  PERFORMS  THE  INITIAL  OBLIQUE  SHOCK  CALCULATIONS  FOR 
0004  C  THE  REGULAR  REFLECTION  SOLUTION.  MN1  IS  THE  STREAMLINE  VELOCITY  IN 
0005  C  FRONT  OF  THE  INCIDENT  SHOCK.  DELANG  IS  THE  FLOWSTREAM  DEFLECTION 

0006  C  ANGLE.  MN2  IS  THE  STREAMLINE  VELOCITY  BEHIND  THE  INCIDENT  SHOCK, 

0007  C  DELMA  ISTHE  MAXIMUM  FLOWSTREAM  DEFLECTION  ANGLE  FOR  REGULAR  REFLECTION. 
0008  C 
0009 
0010 
0011 
0012 
0013 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 


0001 

0002  C- 

0003  C  SUBROUTINE  NONREG  DETERMINES  THE  REFLECTED  OVERPRESSURE  WHEN 

0004  C  REGULAR  REFLECTION  DOES  NOT  OCCUR. 

0005  C - - - - — - - — 

0006 

0007  SUBROUTINE  NONREG(DELHAX,PRESIN,ALPHA1 .MACHN1 ,MACHN2, 

0008  A  PREREF) 

0009  REAL  LOWBND, MACHN1 ,MACHN2,PSIAN0(11) , ANGLIN ( 1 1) ,SLOPES( 11 ) 

0010  REAL  PSISLP(II) 

0011  COMMON/ TRIG/ RTOD.DTOR, PI 

0012  COMMON/REFL/ GAMMA , TOL , UPBND , LOWBND , XI 

0013  COMHON/LAGRAN/ PSI ANG , PSISLP , ANGLIN , SLOPES 

0014  COMMON/  ATMOS/  ATM 

0015 

0016  C - - - 

0017  C  FIRST,  THE  ANGLE  OF  INCIDENCE,  ALPHD1,  WHERE  REGULAR  REFLECTION 
0018  C  STOPS  IS  DETERMINED  BY  CALLING  SUBROUTINE  ALHAX.  NEXT,  THE 

0019  C  ANGLE  OF  INCIDENCE,  STTANG,  WHERE  THE  LINEAR  APPROXIMATION  BEGINS 

0020  C  IS  DETERMINED  BY  SUBROUTINE  LAGRNQ. 

0021  C - — - 

0022 


SUBROUTINE  PRELIM( ANGINC , PRESIN, DELANG, MN1 ,MN2, DELMA) 

REAL  MN1 ,MN2,MN, LOWBND 

COMMON/REFL/ GAMMA , TOL , UPBND , LOWBND , XI 

SINANGsSIN ( ANGINC ) 

MN 1 = ( SQRT ( ( ( GAMMA+ 1 . ) / ( 2 . »GAMHA ) ) »PRESIN+1 . ) ) / 

&  SINANG 
MN=MN1 •SINANG 

DELANGs ATAN ( 1 . / ( TAN ( ANGINC ) * ( ( ( GAMMA+1 . ) * 

&  MN1*®2*.5/(MN**2-1 . ))-1 . ) ) ) 

MN2={ SORT ( ( ( GAMMA- 1 . ) *XI+ ( GAMMA+1 . ) ) / ( 2«GAMMA*XI ) ) ) 

A  /SIN ( ANGINC-DEL ANG ) 

IF  (MN2.LT. 1.)  MN2=1. 

DELMA=(0.7698«(SQRT(MN2#«2-1 .))»*3)/( (GAMMA+1 .)»MN2*«2) 

RETURN 

END 
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0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 


CALL  ALMAX ( ALPHD 1 , PRESIN ) 

CALL  LAGRNG ( ATM* PRESIN , STTANG , PSI ANG , ANGLIN ) 
ALPHD2=ALPHD1-.001 


C - — - - - - - — — - — - - - - - 

C  THE  FOLLOWING  SECTION  OF  CODE  PERTAINS  TO  THE  CASE  WHEN  ALPHA 1 
C  IS  IN  THE  REGION  WHERE  THE  CUBIC  POLYNOMIAL  IS  USED  TO  APPROXI- 
C  MATE  THE  REFLECTED  OVERPRESSURE.  THE  SLOPE,  SL0PE1 ,  AT  ALPHD 1 
C  IS  DETERMINED  BY  FINITE  DIFFERENCE.  THE  SLOPE,  SLOPES,  AT  STTANG 
C  IS  DETERMINED  BY  SUBROUTINE  LAGRNG.  SUBROUTINE  CUBIC  IS  CALLED 
C  TO  DETERMINE  THE  CUBIC  EQUATION  AND  EVALUATE  IT  AT  ALPHA1  TO  OBTAIN 
C  THE  REFLECTED  OVERPRESSURE. 


IF  (ALPHA1. GT.ALPHD2.AND.ALPHA1.LT. STTANG)  THEN 

CALL  PRELIM(ALPHD1 , PRESIN, DELTA, MACHN1 ,MACHN2,DBLMAX) 

ALPHM1 sFALMAX ( MACHN2 , GAMMA) 

CALL  ITRATE ( 0 , ALPHP 1 , DELTA , LOWBND , ALPHM1 ,MACHN2 , PRESIN ) 
QUOTNTs ( 2 . •GAMMA* ( MACHN2»»2« ( SIN< ALPHP 1 ) ) *»2-1 . ) ) / 

&  (GAMMA+1 . )+1 . 

PREFD1 =QU0TNT*XI-1 . 

CALL  PRELIMCALPHD2, PRESIN, DELTA, MACHN1 , MACHN2 , DELMAX ) 
ALPHM2=FALMAX ( MACHN2 , GAMMA ) 

CALL  ITRATE ( 0 , ALPHP2 , DELTA , LOWBND , ALPHM2 , MACHN2 , PRESIN ) 
QU0TNT= ( 2 . •GAMMA* ( MACHN2  »»2» ( SIN ( ALPHP2 ) ) ••2-1 . ) ) / 

&  (GAMMA+1.  )+1, 

PREFD2=QU0TNTfXI-1 . 

SLOPE1s(PREFD1-PREFD2)/(.001«PRESIN) 

CALL  LAGRNG ( ATM*PRESIN , SL0PE2 , PSISLP , SLOPES) 

VALSTTsFLINE ( PRESIN , STTANG ) 

CALL  CUBIC ( ALPHD2 , PREFD2/ PRESIN , SLOPE 1 , STTANG, VALSTT.SL0PE2, 
&  ALPHA 1.PREREF) 

PREREF=PREREF*PRESIN 
END  IF 


C  IF  ALPHA 1 >STTANG,  THEN  THE  LINEAR  APPROXIMATION  IS  USED  TO  DETERMINE 
C  THE  REFLECTED  OVERPRESSURE. 


IF  ( ALPHA1 .QT.ALPHD1 .AND.ALPHA1 ,GE. STTANG)  PREREFa 
4  FLINE( PRESIN , ALPHA 1 ) •PRESIN 
RETURN 
END 


0001 

0002 

0003 

0004 


C  SUBROUTINE  ITRATE  PERFORMS  AN  ITERATIVE  PROCEDURE  TO  DBTBR- 
C  MINE  THE  SHOCK  WAVE  ANGLE  OF  INCIDENCE,  ANGLE,  FOR  A  GIVEN 
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0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

001? 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 


C  FLOW  DEFLECTION  ANGLE,  ANGLE2.  THE  FUNCTION  WHICH  IS  ITER- 
C  ATED  IS  F.  A  SOLUTION  HAS  BEEN  FOUND  WHEN  F=0. 

C - - - - - - - - - - 


SUBROUTINE  ITRATE(L, ANGLE , ANGLE2 , ANGLQW , ANGLHI , 
&  MN2, PRESIN) 

INTEGER  SIGN , SIQNH , SIGNL 

REAL  MN2,ANG(20),VAL(20),L0WBND 

COMMON/TRIG/ RTOD , DTOR , PI 

COMMON/ REEL/ GAMMA , 70L , U  PBND , LOWBNB , XI 

ANGLEL=ANGLOW 

ANGLEHsANGLHI 


C — — — — — — - - - - - — — — — — ^ 

C  THE  RANGE  OF  VALUES  FOR  ANGLE 1  IS  DIVIDED  INTO  20  EQUAL  INTER- 
C  VALS  AND  F  IS  EVALUATED  AT  EACH  INTERVAL.  THE  SMALLEST  VALUE  OF 
C  F  IS  FOUND.  THE  ARGUMENT  ANGLE  FOR  THIS  VALUE  OF  F  BECOMES  ANGLE 
C  AND  THE  ANGLES  TO  THE  LEFT  AND  RIGHT  OF  THIS  POINT  BECOME  ANGLEL 
C  AND  ANGLEH  RESPECTIVELY. 

C - - - - » - - 

ANGINC= ( ANGLEH- ANGLEL ) / 1 9 • 

DO  910  1=1,20 

ANG ( I )  =AN3LEL+  ( I- 1 )  *ANGINC 
V AL ( I ) =F ( L , ANG ( I ), ANGLE2 , MN2 , GAMMA , PRESIN ) 

910  CONTINUE 

VALMIN= 1000000. 

DO  920  1=1,20 

IF  (ABS(VAL(D).LT.VALMIN)  THEN 
VALMIN=ABS(VAL(D) 

IMIN=I 
END  IF 

92Q  CONTINUE 

ILEFT=IMIN-1 

IRIGHT=IMIN+1 

IF  (IMIN.EQ.20)  IRIGHTsIMIN 
IF  (IMIN.EQ.1)  ILEFTalMIN 
ANGLEL= ANG ( ILEFT ) 

ANGLEHa ANG ( IRIGHT ) 

ANGLE={ ANGLEH- ANGLEL)/ 2. +ANGLEL 

C - - - - - - - — 

C  THE  SIGNS  OF  ANGLEL,  ANGLE,  AND  ANGLEH  ARE  COMPARED  TO 
C  DETERMINE  WHETHER  THE  ZERO  OF  F  IS  GREATER  OR  LESS  THAN 
C  ANGLE.  SUPPOSE  IT  IS  DETERMINED  THAT  THE  ROOT  IS  LESS  THAN 
C  ANGLE.  THE  INTERVAL  FROM  ANGLEL  TO  ANGLE  IS  DIVIDED 
C  IN  HALF  AND  ANGLE  BECOMES  ANGLEH  AND  THE  MIDDLE  POINT 
C  BECOMES  ANGLE.  F  IS  EVALUATED  AND  THE  SIGNS  ARE  AGAIN 
C  COMPARED.  THE  PROCESS  REPEATS  UNTIL  ABS[F( ANGLEL )- 
C  F( ANGLEH )]<TOL. 
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0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

00?9 


900  VALUEL  =F(L, ANGLEL , ANGLE2 , MN2 , GAMMA , PRESIN ) 

VALUEHsF ( L , ANGLEH , ANGLE2 , MN2 , GAMMA , PRESIN ) 
VALUE=F ( L , ANGLE , ANGLE2 , MN2 , GAMMA , PRESIN ) 

IF  (VALUEH.LT.O.)  SIGNHs-1 
IF  (VALUEH.GE.O. )  SIGNHal 
IF  ( VALUEL. LT.O.)  SIGNLs-1 
IF  (VALUEL.GE.O.)  SIGNL=1 
IF  (VALUE. LT.O.)  SIGNa-1 
IF  (VALUE. GE.O.)  SIGNsI 
IF  (SIGN.EQ.SIGNL)  THEN 
TEMPsANGLE 

ANGLE= ANGLE+ ( ANGLEH- ANGLE ) / 2 . 
ANGLEL=TS4P 
END  IF 

IF  (SIGN.EQ.SIGNH)  THEN 
TEMP=ANGLE 

ANGLE= ANGLE- ( ANGLE- ANGLEL ) / 2 . 

ANGLEH s TEMP 
END  IF 

DIFsABS ( VALUEH-VALUEL ) 

IF  (DIF.GT.TOL)  GOTO  900 

RETURN 

END 


0001 

0002 

0003 

0004 

0005 

0Q06 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 


C  SUBROUTINE  ALNAX  DTERHINES  THE  MAXIMUM  ANGLE  OF  INCIDENCE,  ALPHD1 , 
C  FOR  ANY  INCIDENT  OVERPRESSURE,  PRESIN. 


SUBROUTINE  ALMAX ( ALPBD 1 , PRESIN ) 

REAL  MACHN1 ,MACHN2 
COHMON/TRIO/ RTOD , DIOR , PI 
DO  1010  1*30,90,5 
ALPHCHspTCR»I 

CALL  PRgLIM( ALPHCH , PRESIN , DELTA , HACHN1 ,HACHN2 ,DELMAX) 
DELMAXsO. 90»DELHAX 
IF  ( DELTA . GT . PSLHAX )  THEH 
DO  1020  JsI-5,1 
ALPHCH =DTOR»J 

CALL  PRELIH( ALPHCH , PRESIN , DELTA, MACH, 1 1 ,MACHN2, 
DF7.MAX) 

DELHAX =0. 90 *DELMAX 
IF  ( DELTA . GT . DELHAX )  THEN 
DO  1030  K=1 ,10 

ALPHCHsDTOR*{ ( J-1 )♦. 1*K) 

CALL  PRELIM  ( ALPHCH, ?RESXNs DELTA, HACHN1 , 
MACH N2, DELHAX) 

DELMAXsO. 9Q*DELMAX 
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0026 

.  IF  ( DELTA , GT . DELMAX )  THEN 

0027 

ALPHD1sDT0R*((  J-1  )+.1*(K-1 )) 

0028 

•  'RETURN 

0029 

•••■■  END  IF 

0030 

1030 

CONTINUE 

0031 

. END  ‘IF 

0032 

1020 

CONTINUE 

0033 

END  IF 

0034 

1010 

CONTINUE 

0035 

END 

0001 

0002  C - - - - - - - - - - - - - - - 

0003  c  SUBROUTINE  CUBIC  FORMS  AND  EVALUATES  THE  THIRD  DEGREE  POLYNOMIAL 

0004  C  WHICH  FITS  THE  MIDDLE  SECTION  OF  THE  REFLECTED  OVERPRESSURE  VS. 

0005  C  ANGLE  OF  INCIDENCE  CURVES. 

0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 


0001 

0002  C - - - - - - - 

0003  C  SUBROUTINE  LAGRNG  PERFORMS  LAGRANGIAN  INTERPOLATION  TO  DETER- 
0004  C  MINE  THE  VALUE,  Y,  OF  A  FUNCTION  AT  ARGUMENT,  X.  THE  ELEMENTS 

0005  C  OF  ARRAY  FUNCT  ARE  THE  TABULATED  FUNCTIONAL  VALUES  FOR  THE  COR. 

000C  C  RESPONDING  ARGUMENTS  IN  ARGUM. 

0007  C - - - - - - - 

0008 

0009  SUBROUTINE  LAGRNG ;X,Y,AFGUM, FUNCT) 

0010  REAL  ARGUM( *  1 ) , FUNCT ( 1 1 ) ,DIFF( 1 1 } ,A( 11} 

0011 


SUBROUTINE  CUBIC(A,D,F,F,I, J,X,Y) 

REAL  I, J,K,L 

COMMON/REFL/GAMMA , TOL , UPBND , LOWBND , XI 
B=A««2/2. 

C=A««3/6. 

G=F»»2/2. 

H=F»*3/6. 

K=G-B 

L=A-F 

COEFF4= ( L* ( I-D+E*L ) +{ J-E) » ( K+A«L ) ) / ( L» ( H~C+B»L ) +K* ( K+A*L ) ) 

COEFF3 = ( ( I-D+E*L ) -( H-C+B»I,  )fiCOEFF4)/( K+AeL ) 

COEFF2=E-A»COEFF3-B»COEFF4 

COEFF1=D-A*COEFF2“B*COEFF3-C*COEFF4 

Y= ( X«*3 )  *C0EFF4/  6 .+ ( X»«2 )  »COEFF3/  2 .  +X»C0EFF2+C0EFF  1 

RETURN 

END 
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FIRST,  THE  FOUR  VALUES  IN  ARRAY  ARGUM  NEAREST  IN  ABSOLUTE  VALUE 
TO  X  ARE  SELECTED  ALONG  WITH  THE  CORRESPONDING  ELEMENTS  OF  FUNCT. 
THESE  ARE  THE  OPTIMUM  POINTS  TO  USE  TO  FORM  A  THIRD  DEGREEE 
LAGRANGIAN  POLYNOMIAL  TO  THE  FUNCTION  AT  X. 


DO  1110  1=1,11 

DIFF ( I ) = APS ( X-ARGUMC I ) ) 

1110  CONTINUE 
SORTEDaO 
LAST= 1 1 

DO  WHILE  (SORTED. NE.1. AND. LAST. GE.2) 

SORTED =1 
LIMIT=LAST-1 
DO  1120  1=1 , LIMIT 

IF  (DIFF(I).GT.DIFF(I+1))  THEN 
SORTEDaO 
TEMP1=DIFF(I) 

TEKP2=ARGUM(I) 

TEMP3=FUNCT(I) 

DIFF ( I ) =DIFF ( 1+1 ) 

ARQUM(I}aARGUM(I+1 ) 

FUNCT ( I ) =FUNCT ( 1+ 1 ) 

DIFF (1+1 )=TEMP1 
ARGUM ( 1+1 )=TEMP2 
FUNCT(I+1)=TEMP3 
END  IF 

1120  CONTINUE 

LAST=LAST-1 
END  DO 
Y=0. 

1125  CONTINUE 

C - - - - - - - - 

C  THE  THIRD  DEGREE  POLYNOMIAL  IS  FORMED  AND  EVALUATED  AT  X  TO 
C  DETERMINE  Y. 


DO  1130  K=1,4 
A  ( K )  =  1 . 

DO  1140  J=1 ,4 

IF  (J.NE.K)  THEN 

A(K) oA(K)*(I-ARGUM( J) )/(AROUM(K)»AROUM(«I)) 
END  IF 
1140  CONTINUE 

Y=Y+A(K)»FUNCT(K) 

1130  CONTINUE 
RETURN 
END 


0001 

0002 

0003 

0004  C - - - - - — — - — - - 

0005  C  FUNCTION  FLINE  FORMS  AND  EVALUATES  THE  LINEAR  FIT  SECTION  OF 
0006  C  THE  REFLECTED  OVERPRESSURE  CURVES.  PRESIN  IS  THE  INCIDENT  OVER- 

0007  C  PRESSURE  AND  ANGINC  IS  THE  INCIDENT  WAVE  ANGLE  AT  WHICH  THE 

0008  C  LINEAR  EQUATION  IS  EVALUATED. 

0009  C - 

0010 

0011  REAL  FUNCTION  FLINE (PRESIN, ANGINC) 

0012  REAL  M,PSIANG( 11) ,PSISLP( 1 1 ) ,ANGLIN( 1 1 ) ,SLOPES( 11) 

0013  COMMON/TRIG/‘RTQD,DTOR,  PI 

0014  'COMMON/LAGRAN/ PSIANG ,  PSISLP ,  ANGLIN ,  SLOPES 

0015  .  .  COMMON/ ATMOS/ ATM 

0016 

0017  C - - - 

0018  C  SUBROUTINE  LAGRNG  IS  CALLED  TO  DETERMINE  THE  SLOPE,  M 

0019  C - - - — 

0020 

0021  CALL  LAGRNG (ATM»PRESIN,M, PSISLP, SLOPES) 

0022  FLINE=M*ANGINC-M«PI/2.+1 . 

0023  RETURN 

0024  END 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 


C - - - - 

C  SUBROUTINE  SCLEAR  ERASES  THE  SCREEN  AND  SENDS  THE  CURSOR  HOME 
C  IF  USING  A  RETRO  VT,  HP2623,  OR  TEKTRONIKS  TERMINAL. 

C - * - - - 


SUBROUTINE  SCLEAR (TRMNUM) 

INTEGER  TRMNUM 
INTEGER* 4  CHANNEL, STAX 
CHARACTER  ESC/ 27/ 

COMMON/IO/CHANNEL 

IF (TRMNUM  ,EQ.  0  .OR.  TRMNUM  .EQ.  4)  THEN 
IAA=LIB$ERASE„PAGE( 1 , 1 ) 
IDAsSCR$SET_CURSOR(fVAL( 1 ) ,*VAL( 1 ) ) 
ELSE  IF ( TRMNUM. EQ. 3)  THEN 
CALL  ASSIGN_CHANNEL 
CALL  OUTPUT(ESC//«H'//ESC//*J») 

CALL  DEASSIGJLCHANNEL 
ELSE  IF  (TRMNUM. EQ. 2)  THEN 
CALL  ASSIGH_CHANNEL 
CALL  OUTPUT (ESC//CHARC 12)) 


75 


0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 


CALL  DEASSIGN_CHANNEL 
DO  1555  1=1,50 

DO  1550  J=1, 10000 
1550  CONTINUE 

1555  CONTINUE 

END  IF 
RETURN 
END 


C - - - — - - - ............... - 

C  THE  FOLLOWING  SUBROUTINES  ARE  USED  TO  CLEAR  THE  SCREEN  ON  THE 
C  HP2623  AND  TEKTRONIKS  TERMINALS  USING  SYSTEM  SERVICE  COMMANDS. 

SUBROUTINE  INPUT (TYPEIN) 

INCLUDE  »($IODEF)» 

CHARACTER* (*)  TYPEIN 
INTEGER* 4  QUAD (2) 

INTEGERS  CHANNEL,  ST  AT,  SYS$QIOW 
DATA  QUAD/2*0/ 

STAT=SYS$QIOW(  ,*VAT.(  CHANNEL)  ,){VAL(IO^TTYRBADALL.OR.IO$HJiOECHO} 
+  , , , ,* REF (TYPEIN) ,*VAL(LBN( TYPEIN) ) , ,QUAD, , , , ) 

IF  (ST AT)  RETURN 

PRINT*,  'READ  ERROR  * ,STAT 
CALL  DEASSIGJLCHANNEL 
STOP 
END 


C - - 

C  SUBROUTINE  OUTPUT 
C - 

SUBROUTINE  OUTPUT(TYPEOUT) 

COMMON/ 10/ CHANNEL 
INCLUDE  ' ($IODEF) * 

CHARACTER* (*)  TYPBOUT 
INTEGERS  CHANNEL, STAT, SYS$QIOW 

STAT*SYS*QI0W( , *VAL ( CHANNEL ) ,SVAL(IO|J#RIT5LBLK.OR.IO$MJJOFORHAT) 
♦  , , , ,*REF(TYPEGUT) ,SVAL(LEN(TYPEOUT) ) , ,*VAL(0) , , , ) 

IF  (STAT)  RETURN 

PRINT*, 'WRITE  ERROR', STAT 
CALL  DEASSIGN_CHANNBL 


0017 

0018 


STOP 


END 


0001 

0002  C — - , — ——— 

0003  C  SUBROUTINE  ASSIGN 

0004  C 
0005 
0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 


0001 

0002  C - 

0003  C  SUBROUTINE  DEASSIGN 

0004  C - 

0005 

0006  SUBROUTINE  DEASSIGNjCHANNEL 

0007  -  COMMON/ 10/ CHANNEL 

0008  INCLUDE  * { $IODEF } * 

0009  INTEGERS  SYS$DASSGN , CHANNEL ,  STAT 

0010 

00 1 1  STATsSYS$DASSGN( %VAL ( CHANNEL ) ) 

0012  IF  (.NOT. STAT)  PRINT*,  'DEVICE  DEASSIGN  ERROR  • ,STAT 

0013  RETURN 

0014  END 


SUBROUTINE  ASSIGN_CHANNEL 
COMMON/ IO/CHANNEL 
INTEGER*4  SXS$ASSIGN , CHANNEL , STAT 
STAT=SYS$ASSIGN( »TT» , CHANNEL, , ) 

IF  (STAT)  RETURN 

PRINT*,  'DEVICE  ASSIGNMENT  ERROR  • ,STAT 
STOP 


0001 

0002  C - - - 

0003  C  FUNCTION  FALMAX  DETERMINES  THE  REFLECTED  SHOCK  WAVE  ANGLE  FOR 

0004  C  MAXIMUM  STREAM  DEFLECTION  BEHIND  THE  REFLECTED  SHOCK. 

0005  C - - - - - - 

0006 

000?  FUNCTION  FALMAX (HACHN2, GAMMA) 

0008  REAL  HACHN2 

0009  SINSQA=( 1 ./(4.*GAKHA*HACHN2**2) )*( (GAMMA+1 . )*HACHN2**2- 

0010  A  3*+0AHHA+SQRT( (GAKHA+1 . )•( (GAHHA+1 . )*MACHN2®*4~ 

0011  4  2 . •( 3 • -GAMMA ) *MACHN2*®2+( GAHMA+9 . ) ) ) ) 
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0012  FALMAX=ASIN ( SQRT ( SINSQA ) ) 

0013  RETURN 

0014  END 


0001 

0002  C - - - - - - - - — — — — - 

0003  C  FUNCTION  F  IS  THE  FUNCTION  EVALUATED  IN  SUBROUTINE  ITRATE  TO 

0004  C  DETERMINE  THE  REFLECTED  SHOCK  WAVE  ANGLE  FOR  A  GIVEN  FLOW  DBFLEC- 

0005  C  TION  ANGLE. 

0006  C - - - - - - - - - 

0007 

0008  FUNCTION  F (L, ARGUM, KNOWN, VEL, GAMMA, PRESS) 

0009  REAL  KNOWN 

0010  VELTEM=VEL 

0011  IF  (L.EQ.1)  VELTEMa(SQRT{((GAMMA+1 . )/(2. •GAMMA) )• 

0012  &  PRESS+1 , ))/SIN(ARGUM) 

0013  F=2 . * ( VELTEM»«2» ( SIN ( ARGUM) ) »*2- 1 . ) / ( TAN ( ARGUM) • 

0014  &  ( VELTEM**2* ( GAMMA+C0S( 2 . * ARGUM ) ) +2 . ) ) -TAN ( KNOWN ) 

0015  RETURN 

0016  END 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 


C  ONESIG  CONVERTS  X  TO  ONE  SIGNIFICANT  FIGURE. 

FUNCTION  ONESIG (X) 

IF  (  X  .GT.  9.  )  THEN 
DO  10  1=1,10 
X*  =  X/10. 

IP  (X  .IT.  9.)  THEN 
X  =  ANINT(X)*(10.eiI) 

ONESIG  a  X 
RETURN 
END  IF 

10  CONTINUE 
END  IF 

IF  (  X  .GE.  .9  .AND.  X  .LE.  9-  )  THEN 
ONESIG  =  ANINT(X) 

RETURN 
END  IF 

IF  (  X  ,LT.  .9  )  THEN 
DO  20  1=1,10 
X  =  X»10. 

IF  (  X  .GT.  .9  )  THEN 
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0026  X  =  ANINT(X)/( 1Q.**I) 

0027  ONESIG  =  X 

0028  RETURN 

0029  END  IF 

0030  20  CONTINUE 

0031  END  IF 

0032 

0033  END 


0001 

C — 

C  TWOSIG  CONVERTS  X  TO  TWO  SIGNIFICANT  FIGURES. 

c— 

0005 

0006 

FUNCTION  TWOSIG (X) 

0007 

0008 

IF  (  X  cGT.  99.0  )  THEN 

0009 

DO  10  1=1,10 

C010 

X  =  X/10. 

0011 

IF  (X  .IT.  99.0)  THEN 

0012 

X  =  ANINT(X)a< 10.*#I) 

0013 

TWOSIG  =  X 

0014 

RETURN 

0015 

END  IF 

0016 

10 

CONTINUE 

0017 

END  IF 

0018 

IF  (  X  .GE.  9-9  .AND.  X  .LE.  99.0  )  THEN 

0019 

TWOSIG  a  ANINT(X) 

0020 

RETURN 

0021 

END  IF 

0022 

IF  (  X  .LT,  9.9  >  THEN 

0023 

DO  20  1=1,10 

0024 

X  a  X»10. 

0025 

IF  (  X  .GT.  9-9  )  THEN 

0026 

X  =  ANINT(X)/(10.**I) 

0027 

TWOSIG  =  X 

0028 

RETURN 

0029 

END  IF 

0030 

20 

CONTINUE 

0031 

END  IF 

0032 

0033 

END 

79 


LIST  OF  SYMBOLS 


am 

a» 

B 

C 

D 

E 

h 

i 

A' 

(Up 

ral 


M 

Ml 

h2 

-*■ 

n 

V 

Pm 

K 

pi 

PR 

T 

+ 

r 


area  of  bore 

propellant  sound  speed  at  muzzle  Immediately  before 
projectile  ejection 

ambient  sound  speed 

fraction  of  propellant  burnt 

propellant  mass 

bore  diameter  of  gun 

internal  energy  of  propellant  gas  immediately  prior  to  projectile 
ejection 

distance  from  the  origin  of  the  contour  plane  to  the  muzzle 
scale  length  for  explosion 

effective  scaling  length  that  varies  with  angle  from  boreline 
projectile  mass 

effective  projectile  mass  accounting  for  bore  friction 
(  «  1.05  mp) 

Mach  number  of  incident  shock 
streamline  Mach  number  In  front  of  incident  shock 
streamline  Mach  number  behind  Incident  shock 
unit  vector  normal  to  contour  plane 

vector  directed  from  the  boreline  to  the  field  point  of  interest 
that  Is  normal  to  the  shockwave  surface 

muzzle  pressure  for  propellant  immediately  before  projectile 
ejection 

ambient  pressure 
pressure  behind  incident  shock 
pressure  behind  reflected  shock 
incident  overpressure  (atm) 
reflected  overpressure  (atm) 

vector  directed  from  muzzle  to  field  point  of  interest 
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LIST  OF  SYMBOLS  (CONTINUED) 
r  magnitude  of  r 

R  gas  constant 

ta  blast  wave  time  of  arrival 

Tm  propellant  temperature  at  muzzle  Immediately  before  uncorking 

Ta  adiabatic  flame  temperature  of  propellant 

Tmean  mean  temperature  of  propellant  gas 
unit  vector  parallel  to  boreline 
U  combined  volume  of  chamber  and  bore 

Vp  exit  velocity  of  projectile 

x  coordinate  having  origin  at  the  muzzle  with  the  direction 

defined  as  perpendicular  to  the  boreline  and  parallel  to 
the  contour  plane.  Shown  in  Figure  1. 

y  coordinate  having  origin  at  the  muzzle  with  the 

direction  parallel  to  the  contour  plane  and  being  In  the  vertical 
plane  encompassing  the  boreline.  Shown  in  Figure  1. 

z  coordinate  having  origin  at  the  muzzle  with  direction 

perpendicular  to  the  contour  plane.  Positive  direction 
corresponds  to  Increasing  distance  from  contour  plane  as  shown 
In  Figure  1. 

Z  (r/t’j1*1 

aj  wave  angle  of  incident  shock 

<*2  wove  angle  of  reflected  shock 

Y  specific  heat  ratio 

6^  flow  deflection  angle  through  Incident  shock 

$2  flow  deflection  angle  through  reflected  shock 

maximum  stream  deflection  consistent  with  regular  shock  reflection 
^  r., 

n  angle  between  r  and  P 

0  polar  angle  from  boreline  to  field  point 

y  momentum  index 

t  vector  along  boreline  designating  apparent  origin  of  blast  wave 
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LIST  OF  SYMBOLS  (CONTINUED) 


ir  cubic  polynomial  used  In  the  approximation  of  the  reflected 

pressure  coefficient  curves 

t  blast  wave  positive  phase  duration 

<j>  angle  between  boreline  and  contour  plane 

x  ratio  of  heat  losses  to  kinetic  energy 

a  roughness  factor 
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1  AFWL/SUL 

Kirtland  AFB,  NN  87117 

1  ASD/XRA  (STINFO) 

Wright- Patterson  AFB,  OH  45433 

1  Air  Force  Armament  Laboratory 

ATTN:  AFATL/DLQDL 
Eglin  AFB,  FL  32542-5000 

1  Commander 

US  Amy  Development  and 
Employment  Agency 
ATTN:  KODE-TED-SAB 
Fort  Lewis,  WA  98433 

1  Director 

Jet  Propulsion  Laboratory 
ATTN:  Tech  Lib 
4800  Oak  Grove  Drive 
Pasadena,  CA  91109 

1  Director 

National  Aeronautics  and 
Space  Administration 
George  C.  Marshall  Space 
Flight  Center 
ATTN:  MS- I,  Lib 
Huntsville,  AL  38512  87 


No.  of 

Copies  Organization 

i  Director 

NASA  Scientific  &  Technical 
Information  Facility 
ATTN:  SAK/DL 
P.O.  Box  8757 
Baltimore/Washington 
International  Airport,  W)  21240 

1  AAI  Corporation 

ATTN:  Dr.  T.  Stastny 
P.0.  Box  126 
Cockeysville,  MD  21030 

1  Advanced  Technology  Labs 

ATTN:  Mr.  J.  Erdos 
Merrick  S  Steward  Avenue 
Westbury,  NY  11590 

1  Aerospace  Corporation 

ATTN:  Or,  G,  Widhopf 
Bldg.  D8  H4/965 
P.0.  Box  92957 
Los  Angeles,  CA  90009 

1  AVCQ  Systems  Division 

ATTN:  Dr,  Q.  Siegelman 
201  Lowell  Street 
Wilmington,  MA  01887 

1  Technical  Director 

Colt  Firearms  Corporation 
ISO  Huy hope  Avenue 
Hartford,  CT  14061 

1  General  Electric  Armament 

8  Electric  Systems 
ATTN:  Mr  R.  Whyte 
Lakeside  Avenue 
Burlington,  VT  05401 

2  Honeywell,  Inc. 

ATTN:  Mail  Station  MN  112190, 
G.  Stilley 
MN  50-2060, 

Mr.  T.  Melander 

600  Second  Street,  North  East 
Hopkins,  MN  55343 


DISTRIBUTION  LIST  (continued) 


No.  of 

Copies  Organization 

1  Hughes  Helicopter  Company 
Bldg.  2,  MST22B 
ATTN:  Mr.  R.  Forker 
Centineila  and  Teale  Streets 
Culver  City,  CA  90230 

1  Martin  Marietta  Aerospace 

ATTN:  Mr.  A.  J.  Culotta 
P.Q.  Box  5837 
Orlando,  FL  32805 

1  AEROJET  Ordnance  Company 

ATTN:  Mr.  A.  Flatau 
2521  Michelle  Drive 
Tustin,  CA  92680 

1  Olin  Corporation 

276  Winchester  Avenue 
New  Haven,  CT  06504 

1  Director 

Sandia  National  Laboratory 
ATTN:  Aerodynamics  Dept 
Org  5620,  R.  Maydew 
Albuquerque,  NM  87115 

1  Guggenheim  Aeronautical  Lab 

California  Institute  of  Tech 
ATTN:  Tech  Lib 
Pasadena,  CA  91104 

I  Franklin  Institute 

ATTN:  Teeh  Lib 
Race  &  20th  Streets 
Philadelphia,  PA  19103 

1  Director 

Applied  Physics  Laboratory 
The  Johns  Hopkins  University 
Johns  Hopkins  Road 
Laurel ,  MO  20707 

1  Massachusetts  Institute  of 

Technol ogy 

Dept  of  Aeronautics  and 
Astronautics 
ATTN:  Tech  Lib 
77  Massachusetts  Avenue 
Cambridge,  MA  02139  8 


No.  of 

Copies  Organization 

1  Ohio  State  University 

Dept  of  Aeronautics  and 
Astronautical  Engineering 
ATTN:  Tech  Lib 
Columbus,  OH  43210 

1  Polytechnic  Institute  of 

New  York  Graduate  Center 
ATTN:  Tech  Lib 
Route  HO 

Farmingdale,  NY  11735 

2  Loral  Corporation 
ATTN:  S.  Schmotolocha 

Sen  Ax ley 

300  N.  Halstead  St. 

P.0.  Box  7101 
Pasadena,  CA  91109 

1  Director 

Forrestal  Research  Center 
Princeton  University 
Princeton,  NJ  08540 

1  Kaman  Tempo 

ATTN:  Mr.  J.  Hindes 
816  State  Street 
P.Q.  Drawer  QQ 
Santa  Barbara,  CA  93102 

1  Southwest  Research  Institute 

ATTN:  Mr.  Peter  S.  Westine 
P.0.  Drawer  28510 
8500  Culebra  Road 
San  Antonio,  TX  78228 

1  Boeing  Company 

ATTN:  C.  R.  Pond 
NS  8C-64 
PO  Box  3999 
Seattle,  WA  98124 

1  FHC  Corporation 

Northern  Ordnance  Division 
ATTN:  Scott  Langlie,  Advanced 
Techniques 

4800  East  River  Road 
Minneapolis  NN  55421 


DISTRIBUTION  LIST  (continued) 


Abe rtieen  Proving  Ground 

Dir,  USAMSAA 
ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen 

Cdr,  USATECOM 
ATTN:  AMSTE-TO-F 

Cdr,  CROC,  AMCCOM 
AT  I'N :  SMCCR-RSP-A 
SMCCR-MU 
SMCCR-SPS-IL 

Dir,  Wpns  Sys  Concepts  Team 
ATTN:  AMSMC-ACW 

Dir,  USAHEL 
ATTN*.  Dr.  Weisz 

Cummings 
Mr.  Garinther 

Dir,  USACSTA 
ATTN:  Mr.  S. Walton 


